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Abstract 

We develop a new 3+1 dimensional Monte Carlo cascade solving the kinetic on-shell Boltzmann 

(N 

^ ' equations for partons including the inelastic gg <-^ ggg pQCD processes. The back reaction channel 

00 ' 

t^^ ■ is treated - for the first time - fully consistently within this scheme. An extended stochastic method 

(N , 

^O i is used to solve the collision integral. The frame dependence and convergency are studied for a fixed 

^^ I tube with thermal initial conditions. The detailed numerical analysis shows that the stochastic 

method is fully covariant and that convergency is achieved more efficiently than within a standard 

geometrical formulation of the collision term, especially for high gluon interaction rates. The 
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(— I ! cascade is then applied to simulate parton evolution and to investigate thermalization of gluons 



X 



for a central Au+Au collision at RHIC energy. For this study the initial conditions are assumed to 
be generated by independent minijets with pT > Po = '^ GeV. With that choice it is demonstrated 
that overall kinetic equilibration is driven mainly by the inelastic processes and is achieved on a 
scale of 1 fm/c. The further evolution of the expanding gluonic matter in the central region then 
shows almost an ideal hydro dynamical behavior. In addition, full chemical equilibration of the 
gluons follows on a longer timescale of about 3 fm/c. 
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I. INTRODUCTION 



The main subject of the heavy ion experiments at the Relativistic Heavy Ion Colhder 
(RHIC) at BNL and at the Large Hadron Colhder (LHC) at CERN is to create a new state of 
matter, the Quark Gluon Plasma (QGP), which is expected to be a transient thermal system 
of interacting quarks and gluons. Due to the confinement free quarks and gluons cannot be 
detected. The search for QGP has to be carried out by analyzing certain proposed hadronic 
and electromagnetic signatures jj, |2, S 14 , S la ID, Isl • However, the possible signatures of 
the QGP may also come in par t from the late time dynamics of a hadron gas formed after 



the phase transition 



layyyyyy 



Therefore one needs detailed informations 
about the creation of the QGP, its lifetime and the hadronization in order to draw reliable 
conclusions. 

Recent measurements jl6| at RHIC of the elliptic flow parameter V2 for semi-central 
collisions suggest that - in comparison to fits based on simple ideal hydrodynamical models 
[ill - the evolving system builds up a sufficiently early pressure and potentially also achieves 
(local) equilibrium. On the other hand, the system in the reaction is at least initially far from 
any (quasi-) equilibrium configuration. To address the crucial question of thermalization of 
gluons and quarks, a number of theoretical analyses have been worked out either using the 



relaxation time approximation 
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Q 
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2l| or performing full 3 + 1 dimensional Monte 



Carlo cascade simulations based on the solution of the Boltzmann equations for quarks and 
,.u„ns H fl Q. TKe «.t parton cascade, VNI, ,^.e, by pQCD i„d„di,, 

binary elastic scatterings (2 ^^ 2) and gluon radiation and fusion (1 ^^ 2) was developed by 



Geiger and Miiller 



22| . In the simulation for a central Au+ Au collision at RHIC energy 



they concluded that a thermalized QGP will be formed at r ^ 1.8 fm/c. However, the onset 
of potential hydrodynamical behavior during the parton evolution was not demonstrated in 
their analyses. In addition, the treatment of the propagation of off-shell partons in their 
approach is not clear from a physical point of view. Recently, Molnar and Gyulassy studied 
he buildup of the elliptic flow at RHIC J2S| applying an on-shell parton cascade, MPC 



2J| (an improved version of ZPC 23/|), in which up to now only elastic gluon interactions 



are included. In their analysis the early pressure can be achieved only if an unrealistic, 
much higher cross section is being employed. Furthermore, it is known that the elastic 
(and forward directed) gg ^^ gg collisions cannot drive the system to kinetic equilibrium, as 



pointed out in Ref. J29|. This would suggest that the collective flow phenomena observed 
at RHIC cannot be described via pQCD. On the other hand, the possible importance of 
the inelastic interactions on overall thermalization was raised in the so-called "bottom up 
thermalization" picture |30|. It is intuitively clear that gluon multiplication should not only 
ead to chemical equilibration j31i], but also should lead to a faster kinetic equilibration 



m 



. This represents one (but not all) important motivation for developing a consistent 
algorithm to handle inelastic processes like gg <-> ggg. 

In solving the transport equations, in most of the cascade models cross sections are 
interpreted geometrically to model the collision processes. It turns out that in dense matter 
when the interaction lengt 
causality violation [2J, 



i Jcr/n is not much smaller than the mean free path of particles, 
will arise in these cascade models and will lead to numerical 
artifacts |36|, |37[. One way to reduce these artifacts is to apply the common test particle 
method (or "particle subdivisions") |38|, |39|, in which the interaction length of the test 



particles is reduced by y/Ntest, while the mean free path is unchanged. Ntest denotes the 
number of the test particles per real particle. However, the limitation of these transport 
models is obvious: Inelastic collision processes with more than two incoming particles cannot 
be straightforwardly implemented since it is in general difficult to determine, for instance, a 
3^2 process geometrically. Therefore, until now, the role of the inelastic processes in the 
formation of the QGP has not been studied fully quantitatively. 

An alternative collision algorithm suggested in 40, |4l|, |4^ dealt with the transition rate 
instead of the geometrical interpretation of cross section and determined proceeding collision 
processes in a stochastic manner by sampling possible transitions in a certain volume and 
time interval. This collision algorithm opens up the possibility to include the inelastic 
collision processes into transport simulations solving the Boltzmann equations 

where C22 and C23 denote the collision term of 2 ^^ 2 and 2 ^^ 3 processes. In this paper we 
will present a newly developed on-shell parton cascade using this sort of stochastic collision 
algorithm. Also the oftenly employed scheme based on the geometrical interpretation of 
cross section is discussed and compared with the stochastic algorithm. In particular, we 
concentrate on the study of the (unphysical) frame dependence. The new transport scheme 
will then be applied to simulate the parton evolution for a central ultrarelativistic heavy 



ion collision at highest RHIC energy. The emphasis is put on the investigation of gluon 
thermalization and their collective dynamics. For this investigation the initial conditions 
are assumed to be generated by independent minijets J43|, |4J]. Other initial conditions, like 
the much discussed "color glass condensate" 45], can also be implemented, but we leave 
this for a future work. For the present study we consider quarks and gluons as classical 
Boltzmann particles throughout the paper. The Pauli blocking and gluon enhancement can, 
in principle, be implemented and will also be discussed elsewhere. 

The paper is organized as follows. In Sec. |n]we consider two-body collision processes and 
contrast the geometrical with the stochastic collision algorithm. The dynamical evolution 
of a system within a fixed box is carried out to study global kinetic equilibration. In 
addition, such calculations are mandatory to debug the operation of the code and to look 
for the limitation of the algorithms. The implementation of the inelastic collision processes 
is described in Sec. IIIII There, we carry out box calculations to study global kinetic and 
chemical equilibration. In Sec. II VI we study thermalization of a parton system in a box with 
initial conditions sampled according to the production of minijets as expected in a central 
heavy ion collision at RHIC. The Lorentz or frame (nondependence) and the convergency of 
results extracted from cascade simulations are investigated in Sec. |3 In Sec. IVII we then 
finally present first results of cascade simulations for a central Au+Au collision at RHIC 
energy. (The readers, who are solely interested in the operation and results of the full 3 + 1 
dimensional cascade, may directly pass to Sec. IVII ) We summarize in Sec. IVIII and give an 
outlook for future works. In Appendices \M and El niore details of the geometrical collision 
algorithm are given. We list the pQCD partonic scattering cross sections in Appendix O 
for two-body processes and in Appendix O for 99 ^^ 999 processes. In Appendix IE] the 
numerical recipes for Monte Carlo samplings are presented. 

II. TWO-BODY COLLISION PROCESSES 

We consider a system consisting of classical, ultrarelativistic particles which are interact- 
ing via two-body collisions. The main emphasis is put on the numerical realization of such 
collision sequences in a relativistic transport simulation, which is theoretically based on the 



solution of the Boltzmann equations (^ with the following collision term given by 

1 f d^P2 1 f (fp'i d% , ,,,, |2/o A4e(4)/ ^ / /N 

-2eJ j2^^f2E2-u I {2.f2E[{2.f2E'J'f'\^''-''''\ ^2^) ^^ Hpi +P2 - 1>, -p,) • 

(2) 

z/ will be set to 2 when considering the double counting if 1' and 2' are identical particles. 
Otherwise v is set to 1. 

Since no mean field is considered throughout the present study, the evolution of particles 
is intuitively straightforward: Particles move along straight line between two collision events. 
After a particular collision the momenta of colliding particles are changed statistically ac- 
cording to the differential cross section. The determination of the collision sequence is, 
however, not unique and depends on the particular numerical implementation. We present 
in this section two numerical methods dealing with the realization of binary collisions. Com- 
parisons between these two methods will be made in detail when investigating kinetic equi- 
libration in a fixed box. We also study any potential (but unphysical) frame dependence of 
transport simulations within both schemes and how to minimize possible deficiencies. These 
results will be presented later in Sec. El 

A. The geometrical method 

In the first method a collision happens when two incoming particles approach as close to 



each other that their closest distance is smaller than a/ 0-22/71", where (T22 denotes the total 
cross section for the colliding particles. In other words, the collision probability is either 
1 or 0, depending on how close the collision partners come together. Since the total cross 
section is interpreted geometrically, we label this procedure the "geometrical method". In 



this picture of the c 
like ZPC bj, MFC 



osest approach,which is already employed in parton cascade models 



2J and PCPC 



collisions do happen one by one as time proceeds. 
The next collision event can be determined by comparing the individual times marking the 
occurrence of the various and possible collisions. 

Unlike the total cross section the closest distance is, however, not invariant under Lorentz 
transformation. This leads to the situation that a particle pair collides in one frame, but 
might not in another frame, which is unphysical. One faces here a violation of covariance. 



which is a historic problem in microscopic simulation within relativistic transport models. 
In the present scheme we define the closest distance in the center of mass frame of the 
individual particle pair [3'J1| and thus make it to be a Lorentz invariant quantity by hand. In 
spite of this definition the covariance of the Boltzmann equation is still not fulfilled, because 



m m- 



the time ordering of collisions might be changed under Lorentz transformation 
Still, for a sufficiently dilute system the geometrical method works rather robust. We will 
continue discussing this problem of covariance violation later in this section and also in Sec. 
IVl Besides the problem just mentioned, the ordering time of one particular collision itself 
which orders the occurrence of all collisions in a particular frame, called lab frame, is not well 
defined. Since we determine the closest distance of two incoming particles in their center of 
mass frame, it is reasonable to define the collision points for the two particles also in this 
frame at the closest distance and at the same time. Consequently both particles, if they do 
collide, change their momenta at the same time in their center of mass frame, but generally 
at different times in the lab frame. (We now denote these individual two times by "collision 
times".) One can now define the ordering time at some stage between these two collision 
times. There is, however, no unambiguous prescription. In general, different choices for the 
ordering time will lead to different collision sequences. This, as numerically verified, does 
not strongly affect the behaviors of physical (ensemble averaged) quantities shown below. 
In our simulation we choose the smaller one of the two collision times as the ordering time. 

n n 

In ZPC [233 and MPC [24^ the ordering time was taken as the average of the two collision 
times. 

In order to demonstrate the correct operation of the numerical realization of the geomet- 
rical method, we will choose a situation when the outcome is known analytically. For this 
purpose we carry out "box calculations" , in which a particle ensemble with a nonequilibrium 
initial condition is enclosed in a fixed box and will evolve dynamically until an appropri- 
ate final time. The collisions of particles against the walls of the box are simply done via 
mechanical reflections. For sufficiently long times, the system should get kinetically equili- 
brated at the end. For a classical, ultrarelativistic ideal gas the energy distribution has the 
Boltzmann form 

^^ ' .e-^/^ (3) 



NE^dE 2T3 
which guides as an analytical reference for the numerical results. The temperature T can 



be obtained from the simple relation between energy and particle density 

e = 3nT , (4) 

where e and n are solely given by the initial conditions. Initially, particles are now distributed 
homogeneously within the box and their momentum distribution is chosen highly anisotropic 
via 

"^^ 6{pT-QGeY)5{p,). (5) 



Ndpxdpz 
In Fig. ^ the final energy distribution from such box calculations for a system of A^ = 2000 

massless particles is depicted. The size of the box is set to be 5 fm x 5 fm x 5 fm. We consider 

isotropic collisions and take a constant total cross section of (722 = 10 mb. The final time is set 

to be 5 fm/c. (As one will shortly realize, this chosen time is sufficient long for the system to 

become equilibrated.) To improve statistics we have collected particles from 50 independent 

reahzations. The dotted line depicted in Fig. [T] denotes the analytical distribution Q with 

temperature T = 2 GeV. We see a nice agreement between the numerical result and the 

analytical distribution except a slight, but characteristic deviation at low energies. We will 

come back to explain this discrepancy immediately. 

Such a successful passing of the previous test is necessary for every collision algorithm, 

but it is still not a sufficient argument to guarantee whether the presented algorithm is 

operating correctly. One has to ask any numerical algorithm for its limitation of correctly 

describing the underlying problem. To be specific when considering the collision integral 

(0), it is not obvious whether the geometrical interpretation of the total cross section is 

a reasonable choice to account for the Boltzmann process. In fact such a description has 

some shortcomings concerning causality violations which have been pointed out for example 

in Ref. J23]- Especially for the algorithm presented above we have to face the fact that 

the collision times of colliding particles are different in the lab frame. This will lead to a 

noticeable reduction of the collision rate compared to one given by the collision integral: 

Assume that the difference of the collision times is Ate. Consequently the particle with 

larger collision time should not collide again during this interval At^, otherwise causality 

would be violated. As pointed out in Appendix for a system in equilibrium the ensemble 

averaged time delay < Ate > depends only on the total cross section and increases with 

the increasing total cross section. This will lead to an artificial increase of the mean free 

path and thus to a decrease of the collision rate. In other words, the collision rate decreases 



when noncausal collisions are forbidden. This problem has also been pointed out in Ref. 
36|,|37|. We can demonstrate this effect employing box calculations, in which we consider 



an initially kinetic equilibrated gas distributed homogenously within the box. The size of 
the box is taken to be the same as in Fig. ^ We employ isotropic collisions with a constant 
cross section of o"22 = 10 mb. In Fig. |21 collision rates are depicted as solid squares for 
several particle densities. The collision rate is obtained here as the time average of the 
collision number. While the box size is fixed, we vary the particle number to get different 
densities. The solid line shows the expected relationship between the collision rate and 
particle density in equilibrium R = no"22 • We see a clear decrease of the collision rate when 



the expected mean free path l/na22 is not much larger than the interaction length W 0-22/71". 
Such a numerical artifact would strongly slow down the kinetic thermalization of an initially 
highly nonequilibrium state, as, for instance, in case of ultrarelativistic heavy ion collisions. 
As also clearly seen from Fig. |2l the collision rate tends to saturate at high density. The 
reason for this is that the collision rate has an upper limit which is exactly the inverse of 
the average collision time difference < Ate > /2 depending only on the total cross section 
as mentioned before. One can compute < Ate > /2 analytically. The detailed calculation 
is given in Appendix 1X1 It turns out that < Ate > /2 = 0.12 fm/c for 0-22 = 10 mb. This 
indicates that the saturation value of the collision rate would be 8.3 fm~^ at high density. 

We now return to the slight discrepancy at low energy as noticed in Fig. ^ and consider 
this as a consequence of the same effect of the relativistic time spread of collisions pointed 
out above, since in this particular situation the particle density is so high that the mean 
free path is one order of magnitude smaller than the interaction length. To confirm this 
suspicion, we carry out similar calculations as in Fig. ^ but with a tiny cross section of 
(J22 = 0.1 mb. The energy distribution, depicted as thick histogram, is shown in Fig. El 
compared with the distribution (thin histogram) obtained by using (J22 = 10 mb. One does 
not see the artificial distortion in the spectrum at low energies any more when the cross 
section and hence the relativistic time spread is small. As a conclusion, the relativistic time 
spread effect not only decreases the collision rate, but also slightly distorts the system out 
of equilibrium. 

To suppress this numerical artifact and hence to conserve Lorentz covariance we employ 



the widely used test particle, or "subdivision" , technique [38|, |39| based on the scaling 

n^nNtest and a ^ a/Ntest, (6) 

where Ntest is the number of test particles belonging to one real particle. While the mean free 



path is unchanged by the scaling, the interaction length is reduced by a factor of ^J Ntest- This 
consequently reduces the relativistic time spread which vanishes in the limit Nt^st -^ oo- The 
open squares in Fig. El denote the results by using Ntest = 50. The tendency of convergency 
towards the ideal limit is visible. 

In Fig. 13 we show the time evolution of the momentum anisotropy defined as the frac- 
tion of the average transverse momentum squared over the average longitudinal momentum 
squared. The initial conditions and parameters are set to be the same as in Fig. ^ The 
dotted line depicts the result without applying the test particle method {Ntest = 1) and the 
dashed line shows the result with Ntest = 50. The results confirm our reasoning that the 
relativistic effect of spreading of the two collision times for a colliding particle pair increases 
the relaxation time for achieving kinetic equilibrium. 

B. The stochastic method 

In the last section we have determined the collision probability of two incoming particles 
by means of the geometrical interpretation of the total cross section. Instead, one can also 



derive t 
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le collision probability directly from the collision term of the Boltzmann equation 



42| . When assuming two particles in a spatial volume element A^a; with momenta 
in the range (pi, pi + A^pi) and (p2, P2 + ^^^2), the collision rate per unit phase space for 
such particle pair can be read off from Eq. (0) 

-Jih 



Atj^A^xA^p, 2E^{2'Kf2E2 



'(2^) 



Expressing distribution functions as 

AN: 

(27r)3 ^ "t^ Pj 

and employing the usual definition of cross section J46| for massless particles 



one obtains the absolute collision probability in a unit box Ax and unit time At 

A A^2^2 A- 

^^^=AiV^AiV^ = "-'^^^A^- ('°) 

Vrei = S/2E1E2 denotes the relative velocity, where s is the invariant mass of the particle 
pair. Unlike in the geometrical method where the collision probability is either or 1, P22 
now can be any number between and 1. (Notice that, in practice, one should choose 
suitable A^x and At to make P22 to be consistently less than 1.) Whether the collision will 
happen or not is sampled stochastically as follows: We compare P22 with a random number 
between and 1. If the random number is less than P22, the collision will occur. Otherwise 
there is no collision between the two particles within the present time step. We call this 
collision algorithm the "stochastic method". Since in the limit At ^ and A^x -^ 
the numerical solutions using the stochastic method converge to the exact solutions of the 



Boltzmann equation 47|, we divide in practice the space into sufficient small spatial cells. 

For a true situation At and A^x have to be taken smaller than the typical scales of spatial 

and temporal inhomogeneities of the particle densities. Only particles from the same cell 

can collide with each other. If a particle pair collides, the collision time will be sampled 

uniformly within the interval {t,t + At). The collision times for both colliding particles are 

here the same. The particle system propagates now from one time step to the next. This is 

different compared to the transport simulation scheme utilizing the geometrical method. 

In general we also might employ, in addition, the test particle technique in order to reduce 

statistical fluctuations of the collision events in cells. Accordingly the collision probability 

is changed to 

o' ^22 At . 

2' '''NtestA^X ^ ' 

by the scaling a -^ a /Ntest- 

In the following we discuss the Lorentz invariance of the stochastic algorithm in the limit 
A^x -^ 0, At ^ and Attest ~^ oo- Since AtA^x, A^p/ AE, the distribution function / 
and the total cross section are Lorentz scalars, it is easy to realize from Eq. ((7j) that the 
collision number AN^^i"^ is a scalar under Lorentz transformations. Furthermore this is also 
true for ANi, the particle number counted within a phase space interval at time t. Hence, 
the collision rate AA^^^^/AAjAr as well as the collision probability P22 are scalars under 
Lorentz transformations. Therefore, in the limit A^x -^ 0, At — i> and Ntest -^ 00 the 
stochastic method yields per se a Lorentz covariant algorithm. However, in practice, a non- 
10 



zero subvoluine A^x and a non-zero tiniestep At disturb full Lorentz invariance explicitely. 
Any potential, but unphysical frame dependence will be discussed later in Sec. El 

To test and demonstrate the stochastic method we again pursue box calculations. The 
initial conditions are the same as in Fig. ^ The size of the box is set as before to be 5 
f m X 5 f m X 5 fm. Since we consider a spatially homogeneous initial situation of particles 
and this configuration will not change very much during particle propagation, we choose a 
straightforward static cell configuration and divide the box into equal cells. The cell length 
is set to be 1 fm in the calculations. We consider isotropic collisions and use a constant total 
cross section of (T22 = 10 mb. Figure El shows the final energy distribution obtained by an 
average over 50 independent runs (with Attest = !)• One clearly recognizes that the stochastic 
collision algorithm also passes this basic test. The agreement between the numerical and 
analytical distribution is perfect and we do not see any distortion in the spectrum in contrast 
to the situation experienced in Fig. |T1 

Since the stochastic method is based directly on the formal collision rate, thus the numer- 
ical realized collision rate should be met in transport simulations if the sampled statistics in 
each cell is sufficiently high. We extract the collision rates from box calculations employing 
the stochastic method and show the results in Fig. IHl as solid squares. The box size and cell 
configuration are set to be the same as in Fig. El The system is taken at thermal equilibrium 
for the initial condition. One nicely recognizes that the squares lay on the expected line. 
(We do mention here that the box size is fixed and we vary the particle number to simulate 
different particle densities. For instance, a density of 1 fm~^ corresponds to a total particle 
number of 125, which means on average one particle per cell. For still lower densities not 
investigated, one would have to work in addition with a suitable amount of test particles.) 

For a system which is initially out of equilibrium the lack of statistics in cells will affect 
the dynamical evolution of the system, since now all cells are correlated during the relax- 
ation time. To study the effect we repeat the same simulations performed for Fig. El starting 
with that particular nonequilibrium initial condition (jH)) and calculate the time evolution of 
momentum anisotropy. We use here the test particle method to control statistical fluctua- 
tions. Figure [71 shows the time evolution of the anisotropy for different test particle numbers 
Ntest- We see that the lack of statistics in cells leads to a slight slowdown in the momentum 
relaxation. This effect is reduced by using larger values for Ntesti which in turn results in 
lower statistical fluctuations. 

11 



Let us summarize with some comparisons between the two simulation methods of treating 
coUisions as presented in this section. In the simulation employing the stochastic method, 
the collision rate is correctly realized if the statistics in the individual cell is sufficiently 
high. In contrast, the collision rate will be numerically suppressed in the simulation using 
the geometrical method, when the mean free path is not much larger than the interaction 
length among test particles. In simulations with both algorithms the test particle technique 
has to be applied in addition in order to solve the Boltzmann equation with sufficient accu- 
racy. For dense and strongly interacting system, convergence of the numerical results with 
increasing test particle number turns out to be more efficient in simulations employing the 
stochastic method than in simulations employing the geometrical method, as shown in Fig. 
0] In transport simulations applying the stochastic method we have to face the difficulty 
of dynamically configurating the space into small cells, which is not necessary in the geo- 
metrical method. Furthermore, the time step has to be chosen much smaller than the cell 
volume to avoid a strong change of the density distribution in cells. This, of course, reduces 
the computing efficiency. In general one should choose such a collision algorithm, so that 
numerical expense is small. However, the stochasic method offers an advanced technique 
when dealing with inelastic collision processes, which is the subject of the next section, 
whereas it might be rather impossible to get a unique and consistent geometrical picture for 
multiparticle transition processes like 2 ^^ 3 for instance. A further comparison between the 
two algorithms will be discussed in Sec. El concerning any potential, but unphysical Lorentz 
frame dependence of the algorithms. 

III. PARTICLE MULTIPLICATION AND ANNIHILATION PROCESSES 

In this section we will now immediately extend the stochastic method to the more compli- 
cated particle multiplication and annihilation processes involving more than two particles. 
These processes are essential to drive the system towards chemical equilibrium and also do 
contribute to kinetic equilibration. The simplest processes are 2 ^-> 3. In physical terms such 
processes will be specified then later in the paper as gluon Bremsstrahlung and its back reac- 
tion. We note that the stochastic method has already been employed for 2 ^^ 3 processes in 
deuteron production pnN ^^ dN J4i| and antibaryon production via, e.g., p + p + co ^^ B + B 
J42J | with much simpler and factorized matrix elements. The true complication in the fol- 
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lowing is to incorporate the true Bremsstrahlung matrix element. Now we will discuss their 
numerical implementations. The implementation of higher order processes is straightforward 
within the extended stochastic algorithm. 

The collision term corresponding the 2 ^^ 3 processes of identical particles is given by 
the expression 

1 1 r d?p2 d^ps 1 f cP'p'i d?p'2 
^^ ^ 2E[2\J (277)32^2 (27r)32E3 217 {2T^f2E[ {2nf2E'2 ^ 

y<f[ n I-Mi'2'^i23|' i2n)U^'\p[+p'2-pi-p2-P3) 

1 r d^p2 1 /■ d^p'^ d^p'2 d^p'^ 

'^2E[J (27r)32E2 317 {2Tiy2E[ (27r)32E^ (2vr)32E^ ^ 

^f[f2f,\Mv2','^12?{2'nY5^'Kp\+p'2+p',-Pl-P2) 

11/" d^p2 d^pz 1 /■ d^p'^ d^p'2 



2Ei 2! J (2Tif2E2 (2t:Y2E^ 2! J {2'nY2E[ {2tt)^2E'2 

x/i /2 /3 |-Mi23^i'2f (27r)^ 5W(pi +P2+P3- P[ - P2) 
1 f d^p2 1 /■ (i^p'^ d^p2 d^p'^ 



2Ei J {2Tif2E2 3! J {2Tif2E[ {2Txf2E'2 {2'nY2E'^ 
x/i /2 |-Mi2^i'2'3f (27r)^(5(^)(pi + p2 - P'l -P'2- P3) ■ (12) 

The collision probability P23 for a particle multiplication process can be derived analogously 

to Eq. (Uni) as 

D ^23 At 



where the total cross section (723 is defined as 
1 1 /■ (i^p'-^ d^p'2 d^p'2 



cr23 



T\Mu^V2'3'\'i2nY5^'\p,+P2-p[~p'2-p'3)- (14) 



2s 3! 7 (27r)32E; (27r)32E; (277)32^3 

One can also extend the geometrical method to the multiplication processes. But it is in 
general impossible to obtain a unified scheme for the annihilation processes in a consistent 
geometrical picture. In contrast, the extension to 3 ^ 2 processes via the stochastic method 
is straightforward. We write the collision rate stemming from Eq. (fT2|) per unit phase space 
in a form like Eq. I^^ 

ANlriVN,,st 1 A^P2 A^P3 /l /2 /3 

Atj^A^xA^Pl 2^1 (27r)32E2 (27r)32E3 Ntest Ntest Attest 

4/(2||^(2|Ie;I^---I'(2^^^^^^^ 

(15) 
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where fi,i = 1,2,3, denote now the phase space density of the test particles. Inserting Eq. 
dHJ into Eq. p5|) gives the colhsion probabihty of a 3 ^ 2 process 






(16) 



for given momenta of the incoming particles in a particular space cell. 1^2 is defined as the 
integral ^ J d^p'^dPp^ ■ ■ ■ in Eg. (fT3j) over the final states. 

Danielewicz and Bertsch J4i| obtained a similar expression for P32 

_ (T12 V3 At 
"^^^-""iV^iV^lA^F' ^ ^ 

when investigating the production of deuterons in a nonrelativistic transport model of low 

energy heavy ion reactions, where they approximately factorized the matrix element into a 

term describing a two-body collision and a term mimicking particle fusion. ai2 is the total 

cross section for the two-body collision and V3 can be interpreted as a volume: Once three 

particles are within this volume, a 3 — » 2 transition may be considered to occur. The volume 

scales with V3 -^ Vs/Ntest when employing test particles. Therefore it is intuitively clear 

why the quantity I32 in Eq. (fTTIjl scales with l/N^^^^. In contrast to Eq. (fTTj) . expression (fTHjl 

is a more general one formulated in a unified manner, and is correct for any given matrix 

elements without any approximations. 

As an example, when considering isotropic 2 <-> 3 collisions for identical particles, integrals 

over momentum space for 0-23 and 132 can be easily calculated analytically and one obtains 



'32 



1927rV23 . (18) 



Applying the probabilities ()13p and (fT^ we are now able to study kinetic and chemical 
equilibration in a box. We assume a system consisting of identical particles and consider 
only isotropic 2 ^-^ 3 collisions. (T23 is set to be 10 mb. As in the box calculations refering 
to Fig. initially the system is chosen to be strongly out of equilibrium according to Eq. 
(0). The particles are distributed homogeneously in the box. The box has a volume of 5 
f m X 5 f m X 5 f m and is divided into equal cells. The cell length is 1 fm. Initially the 
system contains A^o = 2000 massless particles. Newly produced particles will be positioned 
randomly within the individual cells where the transitions occur. Before we come to the 
results, let us determine the final particle density and temperature to be expected when the 
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system becomes thermally equilibrated. For an ultrarelativistic (one component) Maxwell- 
Boltzmann gas the following relations 

rp3 

e = 3neqT and Ueg = — (19) 

hold in equilibrium. One can solve T and Ugg for an energy density given by the initial 
condition. In our case, according to Eq. (jSJ, we obtain T = 1.248 GeV and Ueq = 25.64 
fm~ which is larger than the initial particle density n(to) = 16 fm~ . Figure |H1 depicts 
the time evolution of the particle density obtained from the box calculation. The results 
are obtained by averaging ten independent runs. We see that the particle density increases 
smoothly towards its final value which agrees fully with the analytical expectation. The 
dotted curve presents an estimate made by using the following relaxation approximation 

n{t) = Ueg + {n{to) - rieg) e ^, (20) 

where 6 stands for the relaxation time. In general, for any complex equilibration, this 
quantity will be time dependent. For the estimate the relaxation time is taken by a simple 
fixed value at equilibrium 6 = l/neqCr23 which slightly overestimates the relaxation, as also 
seen in Fig. |H1 In Fig. IHl the final energy distribution is depicted by the histogram. The 
dotted line denotes the analytical distribution with the expected temperature T = 1.248 
GeV. The numerical result agrees again perfectly with the analytical distribution. The fact 
that the final particle density and the final temperature obtained from the inverse slope of the 
energy spectrum are identical to the two analytical values demonstrates that detailed balance 
between the multiplication and annihilation processes is fully realized in our simulations. 
In Fig. El we compare the time evolutions of the normalized particle density (the fugacity) 
and the momentum anisotropy. It turns out that for the given initial conditions the kinetic 
equilibration is slightly slower compared to the chemical equilibration. We notice that the 
quantity 2 <p1> / <p\ >= 2 J d^pp^ f/ J d^pp"^ f is more sensitive to fluctuations than 
n = J d^p/{2TT)^ /, which is the reason why in Fig. El the curve of the fugacity is smoother 
than that of the momentum anisotropy. 

IV. QUARK GLUON PLASMA IN BOX 

A quark gluon plasma (QGP) is suggested as a kineticly and chemically equilibrated 
system of deconfined quarks and gluons. Such state of matter is presumed to have been 
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formed after the big bang and also expected to exist temporarily during the course of an 
ultrarelativistic heavy ion collision in the laboratory. The main goal of the heavy ion collision 
experiments at RHIC and of the future experiments at LHC is to find evidence of such a 
new state of matter, the existence of quark gluon plasma. From the theoretical point of view 
it is also very interesting to address the possibility of the formation of QGP under different 
theoretical assumptions of the initial conditions, and to investigate the further evolution of 
the quark gluon system in space and time. A cascade type transport simulation solving 
relativistic Boltzmann equations for quarks and gluons with Monte Carlo technique is just 
well suited for such a study. Whereas the current parton cascade models, MFC jM], FCFC 
25[ and VNI/BMS 26|, have not included the 2 <-i> 3 processes, we can apply the extended 
stochastic collision algorithm presented in the last section to build up a parton cascade 
describing the space-time evolution of interacting quarks and gluons including gg ^->- ggg 
within the framework of perturbative QCD. As a first application, we restrict ourselves 
in this section to investigate the formation of a quark gluon plasma in a fixed box. The 
convenience is that a thermalized parton system should be formed in any case after some 
time. Although this situation cannot be given in reality, one can still address the way of 
equilibration for different particle species. Furthermore, box calculations offer an essential 
test for the numerical realization of detailed balance of gg ^^ ggg and gg <-> qq processes. A 
realistic space-time approach for the simulation of parton evolution during the early stage 
after an ultrarelativistic heavy ion collision will be presented in Sec. IVII 

The parton interactions include all two-body processes: {1) gg ^^ gg, (2) gg ^-> qq (3) 
gq ^^ gq, (4) qq ^^ qq, (5) qq' <-> qq', (6) qq <-> qq, (7) qq ^-> q'q', and three-body processes (8) 
gg *-^ ggg. The matrix elements squared in leading order of the perturbative QCD are taken 



from Re 
mass 



[Q 



s. |48l . I49l |. We regularize the infrared divergences by using the Debye screening 



2l| rnjj for gluons 

ml = Wnas I T^-{Ncfg + Uff,) (21) 



2 / d^P 1 



and the quark medium mass m^ for quarks 



,2_._. 'v.'-i r i'v 1 



where N^. = "i for SU(3) of QCD and nf is the number of quark flavor. All formulas for 
the differential cross sections are listed in Appendices O and |D| Here we write down only 
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the differential cross sections (or the matrix element squared) of the dominant processes for 
achieving kinetic and chemical equilibration J20|, |31|: 



^^99-99 97ra^ 



dql (gi + m]jf ' 

dq\ 3s(g^ + m?^ ' 



(23) 
(24) 



,2 /'V ' s^ \( 12^^qi 



where g^ = Airas- The matrix element (J25|l describing the (7(7 <-^ (7(7(7 transitions is factorized 
into a part for elastic scattering and a part for gluon radiation (or gluon fusion), q^ and 
kx denote, respectively, the perpendicular component of the momentum transfer and that 
of the momentum of the radiated gluon in the cm. frame. In a dense medium the radiation 
of soft gluons is assumed to be suppressed due to the Landau-Pomeranchuk effect: The 
emission of a soft gluon should be completed before it scatters again. This leads to a lower 
cutoff of kj_ via a step function Q{kj_Ag — coshy), where y is the rapidity of the radiated 
gluon in the cm. frame and A^ denotes the gluon mean free path which is the inverse of 
the gluon collision rate A^ = 1/Rg. Rg is the sum of the rate of the following transitions: 

99 -^ 99, 99 -^ m, 9q ^ 9^, 99 ^ 999, and ggg -^ gg. 

The collision rate is an important quantity governing the time scale of kinetic and chemical 
equilibration. In Fig. ^2 we depict the thermally averged cross section < fj-ezO" > and the 
gluon collision rates as function of temperature for gg -^ gg, gg -^ qq, gq -^ gq, and 
99 ~^ 999 transitions. < Vreic > are calculated numerically, for which we take the screening 
masses obtained at equilibrium {fg = fq = e~^^^) 

m\) = (3 + nA—agT'^ and m^ = -— a^T^ . (26) 

vr 37r 

In the calculations we consider two quark flavors (rij = 2) and employ a constant coupling 
as = 0.3. The corresponding collision rates are obtained hj R = Ug < VreiO' >, where Ug = 
UgT^/Tc'^ is the gluon density in thermal equilibrium. Ug = 2 x 8 denotes the degeneracy of 
gluons. Because of our simple minded inclusion of the Landau-Pomeranchuk effect, the cross 
section cTgg^ggg dcpcuds on the sum of the rates Rg = Rgg^gg + Rgg^gq + Rgg^ggg + Rggg^gg, in 
which, however, Rgg^ggg and Rggg^ggi= Rgg^ggg vo. cquilibrium) depend again on agg-^ggg. 
This problem is solved by a self consistent, iterative computation. Inspecting Fig. ^2 we 
see that the collision rates are proportional to the temperature, which indicates that the 
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< Vrei<y > are inversely proportional to T^. This behavior stems from the fact that the cross 
section (Xgg^gg and cTgg^gq depend mainly on l/rajj and the cross section agg^ggg and CTgg^gq 
mainly on 1/s. Furthermore we realize that the collision rate of the three-body processes is 
in the same order as the rate of two-body gluon collisions. 

We now come to some numerical details when simulating the parton equilibration in a 
fixed box. As shown in Appendix^ the computations of (T23 and 1^2 over momentum space 
are reduced to a four- ()D9|) and a two-dimensional ()D13|) integration respectively. Even 
then, the computations are still time-consuming when cr23 and 132 have to be calculated for 
every gluon doublet and triplet in cells, since the number of integrations is proportional to 
n^ and tt" respectively {n being the total gluon number in an individual cell). In order to 
reduce the computing time, one first thinks of tabulating (723 as well as 132. In simulations 
we then make interpolations using these tabulated data sets. This gives a convenient way for 
obtaining 0-23 because the underlying integral depends on only two parameters, m'jj/s and 
Ag^/s, as mentioned in Appendix O The same data sets have been used for calculating (T23 
in thermal equilibrium as shown in Fig. ^2 In contrast to the case for o"23, 132 depends on 
five parameters (see Appendix EI)- A tabulation of 132 is thus crude due to the limitation of 



herefore we decide to calculate 
with low computing expense 



the storage, which leads to large errors by interpolations. 

J32 in simulations using the Monte Carlo algorithm VEGAS 

(two iterations and 100 function calls). Furthermore, instead of evaluating probabilities of 

all possible collisions, we follow the scheme of Refs. J40, |4l| and choose randomly J\f out of 

the possible doublets or triplets, since in our case the transition probabilities of any channel 

are in fact very small within one time step. In order to achieve the correct collision rate, we 

have to accordingly amplify the corresponding collision probabilities to be 

n{n -l)/2 n{n -l)/2 n{n - l){n - 2)/6 

-^22 -^ -^22 77 , -^23 — > -^23 TT , -f32 "^ -f32 TT ■ ['^ ' ) 

-'V22 ■/V23 ■'V32 

The choices of A/22, A/23 and A/32 are arbitrary. In the following simulations we set A/22 = 

A/'23 = J^32 = n. 

The initial condition for the box calculations is taken by sampling multiple minijet pro- 
duction in heavy ion collisions at RHIC energy ^/s = 200 GeV. Minijets denote on-shell 
partons with transverse momentum being greater than pq, where po is a parameter sep- 
arating the hard, perturbative, from the soft, nonperturbative, nucleon interactions. In 
calculations we set po to be 2 GeV. It had been proposed a long time ago in Ref. J4^ that 
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at RHIC energy the produced minijets take half of the transverse energy. The momentum 
spectrum of the minijets has a power-law behavior and thus the initial condition of the mini- 
jets is strongly out of equilibrium. In the following studies we are interested in the way of 
how thermalization of different parton species proceeds and also interested in the timescales 
of kinetic and chemical equilibration. 

We assume that a nucleus-nucleus collision can be simply modeled as a sequence of binary 
nucleon-nucleon collisions. Then the initial momentum distribution of the produced partons 
is obtained according to the differential jet cross section in nucleon-nucleon collisions J51| 

, 2^^'!j = KY,Xlfa{xupl)x2fb{x2,PT)—^ , (28) 

dp^dyidy2 ^ dt 

where pt is the transverse momentum and yi and 1/2 are the momentum rapidities of the pro- 
duced partons. Xi and X2 are the Feynman variables denoting the longitudinal momentum 
fractions carried by the partons respectively, daab stands for the leading order perturba- 
tive parton-parton cross sections. The phenomenological factor K, set to be 2, accounts 
for higher-order corrections. We employ the Gliick-Reya-Vogt parametrization [52J for the 
parton structure functions faix^p^). For the box calculations we consider gluons stremming 
from a central rapidity region y G [—0.5 : 0.5] as the only initial parton species, since at the 
central rapidity region the partons with small x dominate and these are almost gluons. The 
initial number of gluons is assumed to be 500. 

The primary minijets produced in a real high energy heavy ion collision are distributed 
within a thin disc due to the Lorentz contraction. Instead of such a space-time configuration, 
we assume a homogeneous spatial distribution of partons in the box for simplicity. This 
allows us to still use a static cell configuration. Moreover, all particles are assumed to be 
formed at t=0 fm/c. We will discuss the space-time distribution of the primary minijets 
later in Sec. IVII when considering the parton evolution in a real heavy ion collision. The 
size of the box is set to be 3 fm x 3 fm x 3 fm and the box is divided into equal cells. The 
length of a cell is set to be 1 fm. These settings are tuned as that there will be enough gluons 
(about 15) in each cell during the whole evolution. (For quarks strong statistic fluctuation 
occurs at the beginning of the evolution due to the initial lack of quarks.) 

We employ a constant coupling of a^ = 0.3 in the rest of this section for evaluating the 
screening masses and the cross sections. The screening masses m\, and mi are calculated 
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dynamically according to Eqs. (PT|) and ((221) • The integrations are computed as 



where the sum runs over all particles in a volume V, which should be, in general, small in 
order to maintain the local homogeneity. Since the initial position of partons is distributed 
homogeneously, we extend the sum over all particles in the fixed box. 

The gluon collision rate, which will be employed for evaluating a23 and I32, can be ob- 
tained from the calculated collision probabilities, since the sum of the probabilities of all 
possible collisions gives the average total collision number within the current time step. We 

then have 

y^ pgg^f 
^99'-f = l^\^ ' f = 99^ qq, 999: (30) 

y^ p999^99 
^999^99 = -JJ^^aT ^^^^ 

y^ p9Q^99 
and Rgq^gg = " ^" ^^ , (32) 

where the sums run over possible particle doublets or triplets in the individual cells and also 
over all cells. A''^ denotes the total gluon number in the box. On the other hand, the pP^^^^^ 
and pp99^99 depend again on cr23 and I32 respectively. Therefore, a correct calculation for 
(723 and I32 as well as pP^^^^^ and pP^^^^^ should be a self consistent, iterative computation. 
However, since such computations are too time consuming, we employ the gluon collision 
rate, obtained at the last time step, to calculate (J23 and 132 within the current time step. 

When the parton system becomes fully equilibrated at the later evolution, the final values 
of gluon and quark number should be given by 

iV- = z.,^V^, (33) 

iVf = 2^,^l^, (34) 

where Ug = 2 x 8 and Ug = 2 x 3 x Uf are the degeneracy factors of a gluon and quark 
respectively. The factor 2 in Eq. ()34|) indicates the sum of quark and antiquark. Employing 
the relation 

E = SiN^g" + N^'')T , (35) 

which holds in thermal equilibrium, we obtain the final temperature 
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The total energy E can be determined by the specified initial momentum distribution of 
minijets, Eq. (J2HI)- Considering only up and down quarks (ri/ = 2) we get a final temperature 
of about 430 MeV and thus m^ ^ 0.7GeV^ and m^ ^ O.lGeV^ for a, = 0.3. 

Figure ^1 shows time evolutions of the gluon and quark number. Sixty independent 
realizations are collected to obtain sufficient statistics. We see that the time evolution of the 
gluon number has two stages. At first the gluon number increases rapidly to a maximum 
and then relaxes towards its equilibrium value on a slower scale. The quark number starts 
from zero because of the initial absence of quark species and increases smoothly towards its 
equilibrium value. The gluon and quark number do reach their final values simultaneously. 
These behaviors of Ng{t) and Ng{t) reveal the well-known scenario of two-stage chemical 
equilibration: The gluon system equilibrates at first as if no quarks were there and then cools 
down gradually by producing quark-antiquark pairs until the quarks reach the equilibrium. 
Such two-stage equilibration could also happen in a real high energy heavy ion collision J53| . 

Next we compare the equilibrium values of gluon and quark number of Fig. El with 
the analytical values which one would expect directly from the initial conditions. The final 
temperature in one individual run can be obtained by inserting the total amount of energy 
into expression (j^SI)- Averaged over 60 runs we have < T >= 427.84 MeV. Inserting the 
averaged temperature into Eqs. (fH^ and (jM)) gives < Ng'^ >= 428 and < A^^^ >= 643. The 
values extracted from Fig. ^Jare A^^ = 430 and A^ = 640. We see that the agreements are 
pretty good, which demonstrates that our new cascade algorithm is indeed very successful 
in keeping the detailed balance even for the considered complexity of employing pQCD 
motivated cross sections. We also calculate the equilibrium number of gluons when no 
quarks are considered {nf = 0). In the present situation this is A^'' = 852, which is 
somewhat greater than the maximum of gluon number read off from Fig. ^J since in the 
latter case gluons are already lost due to the production of quark-antiquark pairs starting 
at the beginning of the evolution. 

In Fig. Uni we depict the energy distributions of the partons (gluons and quarks) at 
different times. The initial (t = fm/c) distribution possesses a cutoff at E = po = 2 
GeV and is highly nonthermal. Immediately after the onset of interactions, soft gluons with 
smaller energy do emerge by the process gg — ^ ggg and thermalize very quickly. We see that 
at 0.3 fm/c the energy distribution for partons with smaller energy than 2 GeV is largely 
populated. The hard particles with larger energy are still out of equilibrium. There is still 
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a hump at 2 GeV. This hump will vanish gradually and at 2 fm/c the total distribution 
becomes exponential. One can refer to this stage as the onset of kinetic equilibration. The 
energy distribution at a final time of t = 50 fm/c is also depicted in Fig. ^1 We have 
compared this spectrum to the analytical form Eq. (jH)) with the averaged temperature 
< T >= 427.84 MeV obtained from the initial input. (The analytical distribution is not 
shown in Fig. ^|) The agreement is very good. 

To study the kinetic equilibration in more detail, we calculate the time evolutions of the 
momentum anisotropy 

for gluons and quarks, which are shown in Fig. EI We see that the momentum of the gluons 
and quarks becomes isotropic at almost same time of about 1 — 2 fm/c which is just the 
timescale when the energy spectrum gets exponential, as shown in Fig. ^J However, if one 
looks at the time evolutions of the effective temperatures in Fig. E[ which are defined as 
Tg(t) := Eg(t)/3 Ng(t) and Tg(t) := Eg{t)/3 Ng(t), one notices that between fm/c and 10 
fm/c the temperature of quarks is lower than the one of gluons. The reason is that the quarks 
stem mainly by the gg -^ qq quark pair production and the cross section CTgg^qq is inversely 
proportional to s. Therefore, when the quark production is still more dominant compared 
to the annihilation process, more quark-antiquark pairs with smaller energies are produced 
than those with larger energies, compared to the equilibrated Boltzmann distribution. Cor- 
respondingly, there would be a slight suppression in the energy spectrum of quarks at high 
energy and in the energy spectrum of gluons at low energy during the ongoing chemical 
equilibration. It takes time for the gluon-quark mixture to obtain an identical temperature 
via the gluon-quark interactions. This identical, final temperature is extracted from Fig. 
^Tg = Tq = 429 MeV, and agrees perfectly with the expectation of < T >= 427.84 MeV. 
The parton fugacity is defined as follows 

Xlt):=MM. and AJt) := 3^ , (38) 

where 

Nl\t):=Pg^V and N^^t) := 2Pq^V . (39) 

In Fig. Unithe time evolutions of the fugacity are depicted for gluons (solid curve) and quarks 
(dotted curve). We see that while the gluons approach the chemical equilibrium at about 
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3 fm/c, the quarks do equilibrate later at 20 fm/c. The two-stage chemical equilibration is 
clearly demonstrated in Fig. ^1 

We also depict the time evolutions of the screening masses in Fig. El and of the gluon 
collision rates in Fig. ^J The comparisons of the extracted equilibrium values from the 
figures with the analytical values give perfect agreements. In the small window of Fig. ^J 
the collision rate of gg — > (y'^'^' (upper) and ggg — * gg{\ovfei) are shown by solid lines. We 
see that the two processes occur with the same rates at about 2 ~ 3 fm/c, which is just 
the time scale when the gluons become chemically equilibrated. The identical time scale is 
also obtained from Fig. ^1 We did not depict the time evolution of the rate of ggg -^ gg 
process from 3 fm/c to 50 fm/c, since it is almost identical with that of gg -^ ggg process. 

From the present study of creating QGP in a box some speculations are made when 
we consider parton evolution in a real ultrarelativistic heavy ion collision. (1) Two-stage 
equilibration is a good scenario describing parton thermalization in high energy heavy ion 
collisions. (2) The cross section agg^ggg is in the same order as Cgg^gg and thus the gg ^-> ggg 
processes should play an important role in chemical and as well as kinetic equilibration. 
Analyses based on a full 3 + 1 dimensional transport simulation of the parton evolution after 
a high energy heavy ion collision will be presented in Sec. IVII 



V. TESTING THE FRAME INDEPENDENCE 

The relativistic kinetic equation 

P^dJ = hoii (40) 

is a Lorentz covariant expression. Therefore the covariance of its solution should not be 
affected by the choice of the frame, in which the many-body dynamics is actually described. 
Frame independence must also be fulfilled for any physical observables which can be ex- 
pressed as Lorentz scalars. However, the equation (J40|l cannot be solved exactly in practice 
by applying a transport algorithm. Strictly speaking, the frame independence is not ful- 
filled in any cascade-typ simulations. Our aim in this section is to study potential frame 
dependence in our description employing collision algorithms presented in Sees. O and IIIII 
We will also demonstrate the increasing insensitivity of the particularly chosen frame and 
the convergency of the numerical results when adding more and more test particles into the 
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dynamics. 

As explained in Sec. m the geometrical method is based on the geometrical interpretation 
of the total cross section and the time ordering of the collision events is generally frame 
dependent when the mean free path of particles is in the same order as the mean interaction 
length. In contrast, in simulations employing the stochastic method, which deals with the 
transition rate, a time ordering of the collision sequence is not needed because collision 
events will be sampled stochastically within a time step. Still, one has to be aware that a 
nonzero subvolume of cells and a nonzero timestep disturb the Lorentz invariance. Zh ang 
and Pang had studied already the frame dependence of parton cascade results in Ref. [^ 
applying a parton cascade code with a similar geometrical collision scheme as presented by 
us. They argued that results from parton cascade simulations are not sensitive to the choice 
of the frame when the collision criterion is formulated in the center of mass frame of two 
incoming partons. We will demonstrate the issues in detail in the following considerations 
and calculations. 

A. One dimensional expansion in a tube 

For the purpose of studying the frame dependence we do not need consider a special 
situation. However, as emphasized in the Introduction, the here presented cascade model 
will be applied to simulate the parton evolution in ultrarelativistic heavy ion collisions. 
Therefore it makes sense to consider a one dimensional expanding system as testing ground, 
since at the initial stage of an ultrarelativistic heavy ion collision the partonic system will 
undergo mainly a longitudinal expansion. For convenience, particles of the test system are 
classical Boltzmann particles instead of quarks and gluons. Furthermore, in the present 
section we will employ isotropic collisions and a constant cross section. In order to mimic 
a perfect longitudinal expansion we embed all particles into a cylindrical tube with infinite 
length. The reflections of particles against the tube wall are operated in a same way as 
performed in the box calculations. 



Initially, particles are considered to be therma 
Bj or ken-type boost invariant initial conditions [^ 



in their local spatial element. We use a 



p^ co3h(y — 77) 

/(x,p,r) = e ^M , (41) 
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where r is the proper time r = yW^^z^ and y and t] denote, respectively, momentum and 

space-time rapidity 

1 E + p, I t + z 

?/ = - In — , ri=-m . [Al) 

^2E-p, ''2t-z ^ ^ 

Due to the assumption of the boost invariance, quantities such as particle density n, energy 
density e and temperature T depend only on the proper time r. For an ideal, longitudinal 
and boost-invariant hydrodynamical expansion we obtain 

n(r) = n(ro) ^ , (43) 

T 



4/3 
1/3 



(r) = 6(ro) (^)'', (44) 



nr) = T(ro) (^j . (45) 

Besides the study of the frame dependence we also attempt to address the possibility of 
buildup of an approximately ideal hydrodynamical expansion in cascade simulations when 
the collision rate is considered to be very high. The time dependences (I^Hj) . (HH), and 
fj45|) then serve as ideal references when comparing them with results extracted from the 
numerical simulations. 

To be able to apply the stochastic method, the tube needs to be subdivided into sufficient 
small cells. A static cell structure as configurated in the box calculations is not suitable any 
more for an expanding system. However, since the expansion is only one dimensional, 
we can still employ a static configuration in the transverse plane. Instead of a lattice 
structure, (which will also work,) we make use of the symmetry in the given situation and 
consider a spider web like structure in the transverse plane. Particularly we divide the 
polar angle (p and the radial length squared r^ equally within the interval [0, 27r] and [0, R"^], 
respectively, where R denotes the radius of the cylindrical tube. This division gives a same 
transverse area AF = A</)Ar^/2 for all cells. Longitudinally we have to construct a comoving 
cell configuration which adapts to the expanding system, since, as a reminder, the spatial 
inhomogeneity of particles in the local cells should be small within one time step. Using the 
thermal distribution function (J4H) it can be simply realized by means of the Cooper-Frye 
formula j^ that the particle number per unit space-time rapidity dN/dr] calculated at time t 
in a frame (and also at r as well) is constant, i.e., time independent, when the system expands 
hydrodynamically. This gives us the guideline to divide the tube longitudinally into equal 
small 1] bins. We mark the individual cells [r^j, r^j+i] with the central value r) = {rji + rji^i)/2 
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and the size A?7c = ?7j+i — rji. Then the longitudinal length of a particular cell reads 

Az{t) = t [tanh(r/ + Ar/c/2) - tanh(?7 - Ar/c/2)] (46) 

and increases linearly in time. At time t, when going outwards from the expansion center 
towards the front edges, the cells becomes more and more narrow. Since the particle diffusion 
within a time step should not destroy the homogeneity in the local cells very much, the 
time step has to be chosen smaller than the shortest longitudinal size among all cells. In 
simulations we set the time step to be half of the shortest Az of the cell located at the front 
edge 

At (t) = 0.5 Azminit) = 0.5 1 [t8inh{r]m + Ar]j2) - tanh(r?^ - Ar/e/2)] , (47) 

where rj^ denotes the outermost t] bin. 

With Eq. (j47|) we obtain the collision probability for a two-body process in the central 
cell [t] = 0) 

At 0.5 [tanh(r7„ + Ar/e/2) - tanh(r7„ - Ar/,/2)] 

P22 = Vre^^22^ = Vrel^,, Ax^ 2 tanh(AW2) ' ^'^^ 

For the parameters a22 = 10 mb, Ax±_ = 2.5 fm^, ?]„ = 3.0 and A?]^ = 0.2, the collision 
probability P22 in the central region is expected to be a small value, P22 < 0.004. In order 
to make an estimate of the collision probability in the noncentral cells we go to their local 
comoving frames for convenience, since the collision probability is invariant under Lorentz 
transformations. The time in the local frame of a r^ bin is r = t/7, where 7 = cosh rj denotes 
the Lorentz factor. Suppose that the system undergoes one- dimensional hydrodynamical 
expansion, the collision rate R = n < VreiO'22 > in the local frame of a moving noncentral 
cell is higher than that in the central cell by factor 7, since the particle density is just 7-times 
higher according to Eq. pH|l . [Note that the estimate becomes complicated when the total 
cross section depends on s instead of a constant, since the distribution of s is a function of 
the temperature and the temperatures in the central and noncentral cell are different at time 
t according to Eq. (@S)).] On the other hand the transformed time step At is 7-times smaller 
than At. Therefore the averaged collision number, which is a Lorentz scalar, is the same 
in all cells within a time step At. Furthermore, for the given cell configuration there are 
on average the same number of particles in each cell. This leads to the conclusion that for 
an approximate one dimensional hydrodynamical expansion and choosing a constant cross 
section, the mean collision probability of two incoming particles (for an ensemble average) 
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is the same wherever the colhsion will occur. Due to the fact that the collision probability 
is small we employ the method as explained in Sec. IIVI to reduce the computing time: We 
choose randomly n collision pairs (n being the particle number in a cell) instead of n(n — 1)/2 
possible doublets. The collision probability of each chosen pair is then amplified by a factor 
of [n - l)/2. 

For the numerical simulations we consider a tube with a radius of i? = 5 fm. All particles 
will be produced initially at tq = 0.1 fm/c and are distributed homogeneously within a 
space-time rapidity region 77 G [—3 : 3]. The initial temperature at tq is set to be Tq = 2.6 
GeV and thus the initial particle density is 

^(ro) = 7ri?2 4ro = 1748. (49) 

ar] vr^ 

We have chosen these parameters to achieve initially a dense system. For the cell configu- 
ration we set 

A0 = 27r/8, Ar^ = R'^/i fm^ and Ar/, = 0.2 . (50) 

The transverse area of cells is thus about 2.5 fm and the particle number in one cell is 
around 11. 

The total cross section of the two-body collisions is set to be cr22 = 10 mb if only 2 ^^ 2 
processes are included. We also carry out calculations including both 2 <-i> 2 and 2 <-i> 3 
processes. To be able to make comparisons between simulations without and with inelastic 
processes, we set the cross sections in the latter case to be (T22 = 5 mb and o"23 = 5/2 mb, 
which will lead to the same number of absolute transitions per unit time in both cases. The 
angular distributions of the transitions are considered to be isotropic. 

To study the frame dependence we will simulate the expansion in a so-called lab frame, 
whose origin agrees with the center of the expanding system and in a boosted reference frame, 
which is moving relatively to the lab frame with velocity (3 = — tanhr/o. The situation is 
illustrated in Fig. ^| In the simulations we set t]q = 2. Since particles are initialized 
longitudinally within a limited spatial region in rapidity, the pictures of the expansion in 
the two frames will be quite different. The expansion in the lab frame is symmetric, while 
in the boosted frame the right part of the system expands faster than the left part at late 
times. Therefore the expansion itself is frame dependent at late times due to the limitation 
of the particle initialization. We will concentrate on a so-called central region which is a 
cylinder around 77 = in the lab frame and correspondingly around r/o iii the boosted frame 

27 



with a size of At] = 1. The time evolutions of observables such as ^(t), e(r), T{t), and 
others will be extracted in this central region in the two frames and will be compared. We 
will present the results in Sec. IV CI 

Particles are initialized in the lab frame. At first we sample rj by its uniform distribution 
within f] G [—3 : 3] at the starting time Tq. We then obtain the time and longitudinal 
position of the particle 

^0 = To coshrj, zo = tq sinh rj . (51) 

The transverse positions xq and yo are sampled uniformly within the tube. Finally we 
determine the initial momentum according to the thermal distribution (J4H) at tq for given t]. 
The initial positions and momenta of particles in the boosted frame are obtained by Lorentz 
transformations from the lab frame. 

B. Improved cell configuration 

Before we concentrate on the further analysis, we have to make sure that the cell config- 
uration constructed in the last section is really suitable for an expanding system simulated 
by employing the stochastic method. To demonstrate this we perform a one dimensional 
expansion in the lab frame with the parameters set in the last section and extract dN/drj 
distribution at time t. One expects that the distribution will be constant over a large region, 
since this was the basis motivation for the cell construction. Figure EDI shows the dN/drj 
distribution within an interval of r/ G [—0.3 : 0.3] at time 0.11, 0.13, 0.16 and 0.2 fm/c. 
The dotted line depicts the initial value dN/drjijo) = 1748. Astonishingly, at first sight one 
notices a clear structure in the distribution within the r] bins. (Remember that the size 
of the 7] bins is set to be Arjc = 0.2.) One can also realize that this structure approaches 
a characteristic final shape at late times. The meaning of the structure is that particles 
in a cell are spatially centered. This has no physical reason, but comes from a numerical 
artifact due to the finite size of the cell structure, which can be understood as follows: We 
concentrate on the central 1] bin, t] G [—0.1 : 0.1], and assume that the expanding system is 
in local thermal equilibrium. Any change of the dN/drj distribution in the central r] bin is 
caused from collisions among the particles and from the ongoing particle diffusion. Even in 
the central r] bin the collective motion is still outwards in spite of the small flow velocity. 
There are more particles moving outwards than particles moving towards the center. Sup- 
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pose two extreme cases of collision occurring in the central t] bin: In case 1 two particles are 
moving towards the center and are approaching each other. In case 2 two particles are mov- 
ing outwards and back-to-back. Due to the considered isotropic scattering the momentum 
distribution of the particles after the collision is same in both cases. Since, on average, the 
case 2 happens more frequently than the case 1, one can draw the conclusion that collisions 
in an rj bin tend to bring more particles back into the center than to push them towards 
the outside when the collective flow in an rj bin is indeed directed outwards. This is the 
reason for the artificial structure of the dN/drj distribution in the small rj bins. On the other 
hand, since the distribution of dN/drj is no more constant, the particle diffusion from the 
center outwards is now stronger than the diffusion towards the center. The diffusion is thus 
counterbalancing the particle centralization and the dN/d-q distribution will approach a final 
shape when the balance between the diffusion and the centralization is fully established. 

In Fig. I^we compare the dN/drj distributions at time t = 0.2 fm/c from the simulations 
with At]c = 0.2 and with a smaller size of Arjc = 0.1. In the simulation with Ar]c = 0.1 
we employ 2 test particles per real particle in order to obtain the same statistics as in the 
case with Arjc = 0.2 and Ntest = 1- We see a weaking in the structure of dN/drj, though 
the structure does still exist. In the limit Arjc -^ 0, however, the characteristic substructure 
in the dN/drj distribution will vanish, since the velocity of the intrinsic collective flow in 
the rj bins goes to zero. Therefore decreasing the size of the t] bins and using more test 
particles would be a natural way to reduce this numerical artifact. However, the more test 
particles, the longer will the computing time be. A more elaborate way which does not need 
further test particles is to move the cell configuration randomly from time to time. For 
instance, we move the central rj bin t] G [—0.1 : 0.1] to [—0.1 + ^ : 0.1 + ^], where ^ is a 
random number distributed within [0 : Arjc = 0.2]. Although particles in each t] bin will be 
still centered within each time step after collisions, but because of the random shift of the 
cell configuration the associated center of the bin for a particular particle is also moving, 
so that there is no absolute center for the particle. Therefore, on average, the effect of 
the centralization will be washed out. In Fig. |^we depict the dN/drj distribution from 
simulations employing the improved moving cell configuration with Ar/c = 0.2. We see that 
the distribution is nearly constant and does not show any unwanted substructure. In Fig. 
I7T] we also notice a tiny enhancement of the dN/drj distribution when compared with the 
initial distribution dN/drj = 1748. We will come back to this further artifact in the next 
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subsection. 



Results 



1. 2 <-> 2 processes without test particles 

At first we present the results from simulations without test particles. Figures |221 and 

I2S1 show the time evolution of the particle density, energy density and temperature in the 

central space-time rapidity region in the two frames. The results are extracted from the 

simulations employing the geometrical and stochastic method respectively and are obtained 

by averaging 20 independent realizations. The effective temperature is defined as T = 

e/3n and corresponds to the statistical temperature when the system is at local kinetic 

equilibrium. Otherwise T can be regarded as the mean energy per particle. In the simulation 

with the stochastic method we set the size of the r] bins to be Arjc = 0.2. The time scale 

in Figs. EH and ESI denotes the time in the local frame of the central region. The solid and 

dotted curves depict the results achieved in the lab and boosted frame respectively. The 

thin solid lines show the ideal hydrodynamical limit calculated via a corresponding integral 

of the thermal phase space distribution (J4T|l . Please note that we have taken the size of the 

central region into account. Therefore the hydrodynamical results (J12I), (ISI), and P3|) are 

modified by 

Ar] 
n{t) = aMr = t), a^ := ^ ,,,h(AV2) ^^'^ 

.(«) = a,.(r = t). «.^- 6 tauh(A„/2) /-.,/.*'' (^ + (tanh„T) (cosh,')'" (63) 

T(t) = arT(r = t), 07:=—. (54) 

On 

In the limit Ar/ -^ the additional factors go to 1. 

From Figs. EH and 1221 we see that the frame dependence of the considered quantities is 
quite noticeable in the simulation when employing the geometrical method, while it is rather 
weak in the simulation employing the stochastic method. Moreover and astonishingly, the 
"temperature" in the simulation with the geometrical method is increased at the beginning 
of the expansion. This "reheating" J22| is unphysical, since the isotropic initialization of the 
particle system does not give any reason for an introversive pressure. The gradient of the 
pressure is directed outward, so that in the further evolution the longitudinal work done by 
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the pressure should lead to a cooling of the system. We also rule out any explanation based 
on a possible viscous effect which might bring some effective net energy flow into the local 
region, because there is no reheating in the simulation with the stochastic method. From 
the investigations within a static box we have realized that the collision rate obtained in the 
simulation with the geometrical method will be suppressed when the mean free path is in the 
same order as (or even smaller than) the interaction length. This is indeed the situation at 
the beginning of the expansion in the tube. The suppression of collisions will obviously slow 
down the cooling of the system, but this cannot lead to any reheating. However, the fact 
that particles can interact with each other over a larger distance than the mean free path 
makes it reasonable that the pressure could be incorrectly built up. The effect of the "anti- 
pressure" is thus a numerical artifact. We extract the collision rate and the difference of 
space-time rapidities of colliding particles per collision event < Arj >coU in the central region 
from the simulations carried out in the lab and boosted frame. The results are depicted in 
Figs. 1^ and |2H1 The collision rates are obtained by counting the collision events in the 
central region within a time interval of 0.02 fm/c. It is clearly seen that the collision rates 
in the simulation with the stochastic method agree well with the expectation. The slight 
discrepancy can be understood as the consequence of the relative large size of the t] bins 
(At^c = 0.2). In contrast, the collision rates in the simulation with the geometrical method 
are strongly suppressed at high densities due to the relativistic effect of the time spread of 
the two collision times, as explained in Sec. IIIAl The results of the < Ar/ >coU show that 
particles interact in fact over very large distance at high densities in the simulation when 
employing the geometrical method. The decrease of the < Ar/ >coii at the highest densities 
corresponding to the very beginning of the expansion is due to the fact that at the early 
times particles with large rj are still not formed. In the simulation employing the stochastic 
method the interaction length is, however, controlled by the cell structure. In summary, we 
suspect that the larger interaction distance (compared with the mean free path) may be the 
reason for the "reheating" . 

Figure EEl shows the space-time rapidity distributions at the proper time r = 0.2 and 
1.0 fm/c extracted from the simulations in the lab and boosted frame with the geometrical 
(upper panel) and the stochastic (lower panel) method respectively. The solid (dotted) 
curves depict the distributions in the lab (boosted) frame. The thin solid lines show the 
initial distribution dN/drj^To) = 1748 within r] G [—3 : 3]. We see that the results obtained 
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when employing the geometrical method show a strong frame dependence. A clear hump 
exists around the expansion center r^ = in both frames and broadens gradually. (Note that 
the expansion center in the boosted frame is at r^ = —2 after the shift.) The humps indicate 
a net particle diffusion towards the expansion center, which again can be explained as a 
consequence of the "antipressure" effect: The introversive pressure drives the particles back 
to the expansion center. In the distributions obtained when using the stochastic method we 
see a relative tiny hump at the expansion center which disappears at the later time. The 
slight enhancement has been also noticed in Fig. |^ We recognize that the size of the cell 
bins A?7c = 0.2 is not small enough to overcome the numerical artifact completely. 

In Fig. |27| we depict the momentum rapidity distributions at proper times. The thin 
solid curves show the initial rapidity distribution 

rfiV ^ R^T^Tq sinh(2r7^) 

dy TT cosh.{2rim) + cosh(2?/) ' 

where rjm denotes the boundary of the initial system which has been set to be 3. In the 
upper panel of Fig. |2ZI one also recognizes the particle diffusion towards the expansion 
center, though the effect is not strong. The disributions obtained when using the stochastic 
method show perfect "no frame dependence" and a collective flow outwards to the higher 
rapidity at late times. 

For an initially thermal system it seems reasonable that the system will be still locally 
in or close to kinetic equilibrium during the further expansion. On the other hand, we have 
also realized that numerical artifacts make strong effects at the beginning of the expansion, 
especially in the simulations applying the geometrical algorithm. Therefore it is essential 
to question whether the encountered numerical problem does affect the maintenance of the 
kinetic equilibrium in the cascade simulations of the one dimensional expansion. For this we 
extract the transverse momentum distributions at y = within an interval y G [—0.5 : 0.5] 
at different proper times and compare them with the analytical thermal distributions. In 
Figs. I2H1 and EHl the pt distributions extracted from the simulations in the lab frame are 
depicted. Figure OHl shows the results at r = 1.0 and 4.0 fm/c in the simulations with the 
geometrical method and Fig. 1221 shows the results at r = 0.2, 1.0 and 4.0 fm/c in the 
simulations with the stochastic method. The thermal distributions shown by the solid lines 
are obtained as integral of the thermal particle distribution function (j41|) by means of the 
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Cooper-Frye formula 



1 dN 



N pt dpT Ay 



i'r) = T7 /o \3 A / ^^ / dypTTcosh{y-r])e nr) , (56) 

y=0 ^ (27r)^ AyJ-3 ^-A?//2 

where A?/ = 1 and the temperature T{t) is read off from Fig. |221or Fig. |2niat t = t. We see 
good agreements between the numerical and the analytical distributions, even for the case of 
the geometrical method. The analogous pt distributions, extracted from the simulations in 
the boosted frame (at y = rjo = 2), are also compared with the analytical spectra (both not 
shown in figures). The agreements are perfect as those presented in Figs. EHlandEHl As a 
conclusion, although the expansion does not proceed fully close to ideal hydrodynamics, the 
expanding system is still kinetically equilibrated at least until r = 4 fm/c in the simulations 
with the stochastic method as well as with the geometrical method, although in the latter 
case the cooling of the system occurs much slower. 

As a last point, we show in Fig. EOlthe proper time evolution of the transverse energy 
extracted at y = per unit rapidity from both type of simulations in the lab and boosted 
frame respectively. The thin solid line depicts the result in the hydrodynamical limit 

dEr , , vri?^ /" , ,2 2 ^ 1 \ pt cqsm^-^) 

\J) = -——^dr]dpTPTrcosh{y-r])e ^ij) 



dy 



y=0 



^\R^T'r = \R^nr'J\-'l\ (57) 



The time evolutions of the transverse energy have similar shapes like that of the temperature 
shown in Figs. |221 and |2S1 We also recognize the unphysical "reheating" occurring in the 
simulations with the standard geometrical method. 

Summarizing this section, we have studied the frame dependence of a one dimensional 
expansion in a tube by employing the two collision algorithms presented in this paper. The 
comparisons show that quantities extracted in the simulations with the geometrical method 
have a much pronounced and unphysical frame dependence. Numerical artifacts are very 
significant in these simulations, especially at the beginning of the expansion when the system 
is very dense. In contrast, the results obtained from the simulations when employing the 
stochastic method show almost "no frame dependence". 
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2. 2 <-> 2 'processes with test particles 

The time evolutions of the particle density, energy density and temperature depicted in 
Figs. 1221 and 1221 demonstrate that simulated dynamics does not undergo an ideal hydrody- 
namical expansion. On the one hand, it is true that the ideal hydrodynamics cannot be 
realized in simulations with finite collision rate. One has to take the finite viscosity into 
account. Thus it is interesting to make comparisons between the transport results and those 



calculated from viscous hydrodynamics 
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58|,|59|. This subject is, however, beyond the 



scope of this paper. On the other hand, even the viscous expansion cannot be solved ex- 
actly due to the limitation of the numerical implementations. Especially, as observed in the 
simulations with the geometrical method, the numerical artifacts make strong unphysical 
effects. In this section we introduce the test particle method into the dynamics to reduce 
this numerical uncertainty and to study the convergency of the transport solutions. 

From the experiences in the box calculations, one realizes that the computing becomes 
more time consuming when more and more test particles are added into the simulations. 
One way to reduce the computing time in the present case is to consider a tube with a 
smaller radius. The (real) particle density is however unchanged. In simulations with the 
geometrical method we set the radius of the tube to be R' = R/y^Nfest with R = 5 fm. 
However, in simulations with the stochastic method we instead keep the radius of 5 fm, in 
order to be able to refine the cell configuration. 

Figure EH depicts the relative frame dependence of the particle density, energy density, and 
temperature extracted in the central region in the simulations with the geometrical method 
with Ntest = 1, 4, and 25, respectively. The simulations are performed in the lab frame. We 
obtain the results by averaging 20, 2, and 20 independent realizations, respectively. Note 
that the simulation with Ntest = 4 is exceptionally carried out with the default radius of 
R = 5 fm. We see that the potential frame dependence is more and more reduced when 
more and more test particles are considered. The reduction of the frame dependence is also 
clearly demonstrated in Fig. ESI Here the distribution of the space-time rapidity obtained 
with Ntest = 25 is compared with the distribution without test particles (or Ntest = 1) at 
r = 0.2 fm/c. The hump, which exists in the distribution without test particles due to 
the artificial back diffusion, does not occur with Ntest = 25. For the case employing the 
stochastic method it is not necessary to study the reduction of the frame dependence with 
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the test particle method, since the frame dependence is actually very weak even without 
test particles (see Fig. 1^ . 

We also employ the test particle method to study the convergency of the transport solu- 
tions. Figure inSl shows the time evolution of the temperature extracted in the central region 
in the simulations with increasing test particles in the lab frame. The size of the r] bins 
constructed in the simulations with the stochastic method is refined to Arjc = 0.2/Ntest- 
There are on average 11 test particles in one cell. (We have also performed simulations 
with doubled test particle number in one cell to increase the statistics. The outcome shows 
almost no changes.) From Fig. ESI we see the clear tendency of convergency. The time evolu- 
tion of the temperatures extracted from the simulations with the geometrical and stochastic 
method converge towards almost the same curve. However, it is obvious that the solution 
obtained with the stochastic method converges more efficiently than the solution obtained 
with the geometrical method. Therefore, we do favour the stochastic method to be applied 
in transport simulations of system with high particle density. Furthermore, we see that the 
effect of the artificial reheating, appearing in the simulation with the geometrical method 
with Ntest = 1, reduces and vanishes in the simulations when employing higher test particles. 

In Fig. OUwe depict the collision rate and the mean difference of the space-time rapidities 
of colliding particles per collision < At] >coii in the simulations with the geometrical method 
with increasing test particles. We see that the collision rate increases when using more 
test particles. However, even for Ntest = 900 the collision rates at high densities are still 



suppressed. The reason is that the interaction length decreases only with Ij y/ Ntest- We 
also see that the < Ar/ >coii decreases when the number of the test particles increases. 
Putting Fig. 1^ in relation to Fig. 1221 confirms our suspicion in the last subsection that 
unwanted collisions over large distances may lead to the buildup of "antipressure" which 
then influences the particle diffusion. We mention that the same numerical artifact has been 
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found in the studies in Refs. 



|2j,|37| 



3. Including 2 •s-^ 3 processes 

We now include the inelastic 2 ^^ 3 processes into the dynamics of the one dimensional 
expansion in the tube and study the frame dependence for the new situation. The stochastic 
method is applied to simulate the (in)elastic collisions whose cross sections are set to be 
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(722 = 5 mb and a23 = 2.5 mb. These parameters lead to the same rate of the elastic and 
inelastic transitions. We consider isotropic collisions and set the size of the rj bins to be 
Ar], = 0.2/iVte.t. 

In Fig. 123 we show the time evolutions of the particle density, energy density, and 
temperature extracted in the central space-time rapidity region from the simulations with 
Ntest = 1 carried out in the lab and boosted frame. The results are absolutely frame 
independent. Comparing to the results with only two-body collisions shown in Fig. |22l 
we notice that the particle density is slightly enhanced. This enhancement is not due to 
any numerical artifacts, but the consequence of the chemical equilibration: In the thermal 
equilibrium the particle density is related with the temperature by n = T^/tt^. Since 
during the expansion the temperature is always higher compared to the ideal hydrodynamical 
limit due to finite viscosity, therefore, there have to be more particles being produced than 
annihilated in order to account for the undersaturated system and to achieve a new balance. 
To address the chemical equilibration we concentrate on the time evolution of the fugacity 
which is defined as X{t) = n{t)/rf'^{t), where 

n^«(t) = a„ n^'^{T) = a„ ^ = ^ ^ • (58) 

an and a-p are factors given in Eqs. (|52j) and (j54|) taking the size of the central region into 
account. The T(t) in Eq. (jHH|) is just the extracted temperature from the simulation. Figure 
Oni depicts the time evolution of the fugacity. We see that the chemical equilibrium is almost 
achieved and maintained during the expansion in both frames. We have also extracted the 
Pt distributions and compared with the analytical spectra at different times. The results 
show that the kinetic equilibrium is also maintained during the expansion. 

The collision rates of 2 ^^ 2, 2 ^ 3, and 3^2 processes, extracted from the simulation 
in the lab frame, are depicted in Fig. IHZI We see perfect agreements of the extracted 
collision rates with the expectations. Furthermore, the collision rates of 2 — i> 3 and 3 — >■ 
2 processes are almost identical, which demonstrates once more the maintenance of the 
chemical equilibrium in the expanding system. 

We show in Fig. IHHlthe particle distributions versus the space-time rapidity rj and versus 
the momentum rapidity y at t = 0.2 and 1.0 fni/c, obtained from the simulations with 
Ntest = 1 in the lab and boosted frame. The frame dependence is not noticeable and lies 
within the statistical errors. In Fig. EHlwe also see the enhancement in the particle number 
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over a large range due to the slight particle production in the ongoing chemical equilibration. 

Finally, though the results from the simulations above are largely frame independent, 
convergence to the correct Boltzmann solutions requires using more test particles. Due to 
the settings of (722 = 5 mb and 02^ = 2.5 mb one would expect that the total collision rate 
including elastic and inelastic processes is the same as that in the simulation with purely 
elastic collisions and cr22 = 10 mb. Therefore, the convergence with increasing test particles 
would be exactly the same in both cases. Still, as realized from the above comparison 
between Figs. EH and 123 the temperatures (and the number densities as well) in the central 
space-time rapidity region are slightly different due to the new balancing as explained above. 
Therefore, the convergence of the temperature, for instance, will not be the same as that 
shown in the lower panel of Fig. 1221 On the other hand, the new balancing should not 
affect the time evolution of the energy density. (This quantity is shown in Fig. 1221 for the 
pure 2 ^^ 2 and the 2 ^-i> 2 + 2 ^^ 3 case with increasing Ntest-) We see the exactly same 
convergence of the energy density in the central region. 

After this exhaustive discussion of testing the operation of the cascade, we now proceed 
to describe real heavy ion collisions. 

VI. FULL 3 + 1 DIMENSIONAL OPERATION OF THE PARTON CASCADE FOR 
CENTRAL AU+AU COLLISIONS AT RHIC: KINETIC AND CHEMICAL EQUI- 
LIBRATION 

In this section we take the step to simulate the space time evolution of partons produced 
in a central Au+Au collision at maximal RHIC energy i/i = 200 GeV by means of the 
well-tested stochastic collision algorithm. The simulation is performed in the center-of-mass 
frame of the colliding nuclei. For the present and first exploratory study we include only the 
pQCD motivated gluonic interactions gg ^^ gg and gg ^^ ggg in the dynamical evolution. 
Simulations with all parton degrees of freedom will be postponed to a sequent paper. 

A. The initial conditions 

The initial conditions for the parton cascade are assumed to be generated by minijet 
production in a central Au+Au collision modeled via multiple, binary nucleon-nucleon col- 
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lisions |43|, |4J|. Of course, this is a strong assumption. The picture of the very early stage 
of the colhsion, when potentialy the partons are freed from the two nuclear wave functions 
and do become on-shell particles, is crucial for all kinetic cascades which can only describe 
the further evolution. Hence, one has to incorporate a physical model for describing the 
very initial phase of liberated partons serving as initial condition for cascades. One such 
physical picture is based on the idea of a free superposition of minijets being produced in 
the individual semihard nucleon-nucleon interactions. Another and much celebrated picture 

n 

is the so-called McLerran-Venugopalan model or color glass condensate |45|| , which is based 
on the idea of gluon saturation of the QCD structure function of the nuclei at sufficiently 
low X. The so-called bottom up scenario 30j of thermalization relies on these initial con- 
ditions, where in the leading order of the coupling constant a^ the various time scales of 
kinetic evolution is parametrically estimated. We will leave this as an important task for fu- 
ture investigation. At present, we choose the liberation of minijets as the initial conditions. 
Minijets denote the on-shell partons with transverse momentum being larger than a centain 
cutoff value of po ~ 2 GeV. Since no nuclear effects like shadowing at small x are considered 
in the present study, the averaged number of produced partons is then just proportional to 
the number of binary nucleon-nucleon collisions 

< N.parton > = 2 a,et Taa(6 = 0) . (59) 

TAA{i> = 0) denotes the nuclear overlap function for a central nucleus- nucleus collision. The 
overlap function is given by 

TabC^) = / d^XTidzi d^XT2dz2 nA(ri) ^^(ra) (5^(b - (xn - x^a)) • (60) 

n^/B(r) is the nuclear density. In physical terms, aTAsip), where a denotes the total nucleon- 
nucleon cross section, gives roughly the number of binary semihard nucleon-nucleon collisions 
in a. A + B collision at impact parameter h |44:]. The total jet cross section aj^t in a nucleon- 
nucleon collision at ^/s = 200 GeV is calculated by integrating the differential jet cross 
section (J2HI)- The factor 2 in Eq. (jSHl) indicates that minijets are produced in pair. Employing 
the Woods-Saxon distribution for the nuclear density of a Lorentz contracted nucleus 

7 no 
UAi^^Tl, Zi) = — , ^ , (61) 



l + Exp[{^xl, + {-fz,y-RA)/d 

where d = 0.54 fm, Ra = 1.12A^^^ — 0.86A^^^^ fm, and uq is determined from the normaliza- 
tion / d^ri ua = A, one estimates that with a cutoff po = 2 GeV about 1200 minijets will be 
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produced in a central Au+Au collision at maximal RHIC energy, in which 70% are gluons. 
We note that this number does crucially depend on the cutoff po, which makes the minijet 
picture not so much promising. On the other hand, one might improve this by choosing 
some self-consistent relation for this crucial parameter J60|. 

The initializations of the individual produced partons in space-time and in momentum 
space are then realized statistically as follows: The momenta are sampled according to the 
differential jet cross section (j^Hj) . This sampling has already been performed in Sec. IIVI 
when the thermalization of a parton system was studied within a fixed box. The space- 
time coordinates of the partons are obtained by a simple geometrical picture when the two 
Lorentz contracted nuclei do overlap. For convenience for the moment, we set the zero point 
of the time scale to be the moment of the full overlap. Then the longitudinal positions of 
the two nucleus centers are then at ±vt, respectively, where v is the velocity of the nuclei. 
One now identifies the intrinsic coordinates zi and Z2 in Eq. ()60p with the global space and 
time coordinate 

zi = z — vt and Z2 = z + vt . (62) 

Changing from Zi and Z2 to z and t yields for Eq. (|6(J|) for 6 = and A = B 

TaaO^ = 0) = cfxTid'^XT2'^vdtdznA{^Ti,z — vt)nA{'x.T2,z + vt)5^{xTi—x.T2) 

= d'^XTi2vdtdz nA{x.Ti, z — vt)nA{'x.Ti, z + vt) . (63) 

One thus receives the statistical distribution for sampling the space-time coordinates of the 
individual produced partons 

— - — ^"'^*°" nAixTi,z-vt)nA{yiTi,z + vt). (64) 

d^XTi dz at 

The probability for producing a parton at (r, t) is thus proportional to the convolution of 
the nuclear densities of the two overlapping nuclei at the individual space-time point. Due 
to the choice of the zero point in time, about half of the produced partons are liberated 
at negative times. Therefore, with this convention of the zero point in time the space-time 
rapidity 7] (j^^ is not a well-defined quantity. In order to correct this, we shift all the 
times to be larger than the absolute values of the corresponding longitudinal positions, i.e., 
t — > t + tg > |z|, with a uniquely chosen t^. This actually implies that t^ ~ Q.^Ra/Ii i-e., 
half of the overlapping time. Since we apply Woods-Saxon distribution (jM|) for the nuclear 
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density, we cannot exactly specify when the first touch of the two colhding nuclei occurs, ts 
is thus - strictly speaking - a parameter in our simulation. For a larger tg particles pile up in 
the central space-time rapidity region and for a smaller ts particles distribute within a wider 
rapidity range during the very early evolution. On the other hand, independently on the 
chosen ts, most of the partons are in fact produced in the central rapidity region due to the 
geometry of the overlapping nuclei. In the simulations we determine ts with the assumption 
that the initial partons are distributed within a space-time rapidity range of r/ G [—5 : 5]. 

In the above picture concerning the implementation of the space-time production of mini- 
jets, it is assumed that partons become immediately on-shell when the (semi)hard nucleon- 
nucleon collisions occur. Alternatively, one may introduce an additional formation time for 
every minijet, Atf = coshyArf ^ coshy ■ l/pr, which models the off-shell propagation of 
the freed partons. Within that time span, one assumes, for simplicity, that the still virtual 
parton does not interact and moves with speed of light. We realize, confirmed by numer- 
ical simulations, that the introduction of such a formation time does not affect our main 
findings too much. A further brief discussion will follow in Sec. IVIDI We note that the 
consideration of the initial conditions of partons is a reasonable description of the minijet 
production in space-time according to the overlapping of two heav y io ns. These are different 
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28|. Therefore, the initial 



from the Bjorken-type initial conditions as considered in Refs. 

correlation between the space-time and momentum rapidity, f] — y correlation, cannot be 

expected as the simple rj = y. Detailed analysis will be shown in Sec. IVIDI 

B. Cell configuration 

To be able to apply the stochastic method to simulate the full collision sequences, we 
divide the space into appropriate cells. The individual cell structure has to be considered 
selfconsistently to be well suited to the details of the dynamical evolution of the parton 
system. Since it is not clear whether the parton evolution is invariant under the Bjorken 
boost, a configuration with constant division in space-time rapidity. At] =constant, as chosen 
in Sec. |3 when simulating one dimensional expansion of a thermal system in a tube, is 
here not really reliable. Cell structure should be refreshed every time step to adapt to 
the dynamical parton evolution. In principle, one dimensional expansion is still a good 
approximation for the whole parton evolution for the first few fm/c after a nucleus-nucleus 
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collision. We thus still employ a static cell configuration in the transverse plan: Cells are 
transversely set as squares with a length of 0.5 fm. Longitudinally, space is divided into Az 
bins, where each bin contains about the same number of test particles. This ensures the 
same statistics for each bin and automatically adapts to the density profile of the evolving 
parton system. This dynamical structuring begins at the center of the fire-ball and then 
proceeds to the very outside. Test particles from the far outside are not included into the cell 
configuration, because there the density distribution is too inhomogeneous. Instead, we then 
consider only elastic scatterings among these partons treated via the geometrical method. 
To obtain sufficient statistics, one has to tune the test particle number in each bin to be large 
enough: It turns out that a number of 20 test particles on average in each cell is sufficient 
during first 4 ~ 5 fm/c. However, in the region with lower particle density, especially in the 
transversal surface, there are not enough test particles. If the test particle number in a cell 
is less than a certain cutoff, which is set to be 4 in the simulations, we treat test particles 
in this cell again only by means of elastic scatterings with the geometrical method. How 
fine the longitudinal bins would be, depends on Ntest, the number of test particles per real 
particle. We set Ntest = 60 in the simulations. In total, this leads to an equivalent division 
of roughly equally sized bins in space-time rapidity with Arj = 0.1 ~ 0.2. Furthermore, in 
order to avoid that particles belong to the same cell for too long time, as discussed in Sec. 
fVl we randomly shift the cell structure by a small amount in the longitudinal as well as 
transversal direction after every time step. 

Besides this fine mesh of cells we also have to choose a sufficiently small time step to 
prevent a too strong change of the spatial configuration in each local cell. In the simula- 
tions, this time step is time dependent and is determined to be the one-fifth of the smallest 
occurring cell length. For the case that a collision probability turns out to be greater than 
1, all operations done within the current time step are redone with an appropriately chosen 
smaller time step. 
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C. Assumptions 

We calculate the dynamical screening mass m|, in a similar way as done for the box 
calculations in Sec. |TV| [20| 



771, 



ij^ - 167ra, / Nc fg ^ IGirasNc— 2^ — . (65) 

The evaluation is carried out (quasi)locally. V denotes the volume of a local region and the 
sum runs over all test particles in the region. The presence of the cell structure makes it 
reasonable to calculate the screening mass in each cell. However, the statistical uncertainty 
due to fluctuations is still large, since there are at maximum 20 ~ 30 test particles in one 
individual cell, and thus an extraction of the particle phase space density / is not precise. 
If one assumes that the expansion in the first few fm/c is mainly longitudinal, and further, 
that the transverse parton distribution is homogeneous over a large transversal area, one 
can extend the sum in Eq. (pTK|) over a more broader region compared to the individual 
cell. In the simulations we consider a volume ^ as a cylinder with a radius of 6 fm in 
the individual Az bin. Within each bin m]^/as is assumed to be transversely constant. 
This approximation will lead to an underestimate of m'jj/as in the very central area and 
an overestimate in the outside area when the transverse flow builds up, since within the 
same Az bin the particles moving with larger transverse velocity have larger energy and 
thus make a smaller contribution to the sum in Eq. (jU3j) than the particles moving with 
smaller transverse velocity. The radius is a parameter which we set to be roughly equal to 
the radius of a Au nucleus. It turns out that the influence of this parameter on the screening 
mass is still quite sensitive at least for late times. A future improvement will be to simulate 
the parton evolution within a parallel ensemble technique, which will give the possibility to 
extract the particle phase space density locally more precise. 
The coupling a^ is assumed to be 

as(<5^) = as{s) = 7^^ — - — . , . ,.2 r (66) 

for individual collisions, where s denotes the invariant mass of a particular colliding system 
of two or three particles. We set the quark flavour Uf to be 3 and Aqcd to be 200 MeV. 
In general, Q^ in Eq. (|66|) stands for the momentum transfer in collision such as in deep 
inelastic scattering. For many-body collisions, however, the scale Q^ is not unambiguous. 
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The gluon collision rate Rg, which will be employed to determine the effectively in- 
corporated Landau-Pomeranchuk suppression in the gg ^^ ggg processes by means of a 
low-momentum cutoff, is evaluated locally in cells 

^a ~ ^gg^gg "•" ^gg^ggg "•" ^ggg^gg ' 

where 

y^. pgg^f 
^99^^ = In\t ' f = 99,999, and (67) 

y^ ppgg^gg 
Rggg^gg = \^^^^ ■ (68) 

The PjS denote the respective individual collision probabilities. The sum over Pi gives 
the mean number of collisions occurring during a time step At in a cell with A^^ gluons. 
At denotes the corresponding time interval in the comoving frame Ar = At/7, where 
7 = l/\/l — v^/c^ and v is the collective velocity of the moving cell. For a cell with about 
20 gluons there are totally 200 individual possible gg —>■ gg and gg —>■ ggg collisions each 
and 1200 possible ggg —>■ gg collisions. The statistics is high enough to ensure evaluations 
of the collision rates in local cells to be sufficiently precise, in contrast to the calculation of 
the screening mass. Still a problem remains, which is that the calculation of the collision 
probability for a ggg -^ gg process by a two dimensional integral is time consuming. To 
reduce the computing time we have to take the following approximation, which has already 
been applied for the box calculations in Sec. IIVI We randomly choose about 20 gluon triplets 
instead of the total 1200 combinations and compute the amplified collision probabilities ()27j) . 
Therefore the statistical fluctuation of the collision rate Rggg^gg is stronger than that of the 
others. Also, when extracting the velocity of an individual space element we encounter the 
same difficulty of insufficient statistics as explained by calculating the screening mass. We 
assume that all the cells in a Az bin have the same longitudinal velocity 

V f '^^P 2i f 1 

^/(0/ Ng^E^ 

where the sum runs over the test particles within a cylinder with a radius of 6 fm in the 
considered Az bin, and A'"^ denotes the gluon number in the cylinder. The transverse 
component of the velocity is set to be zero. In principle, this assumption can be corrected 
when a parallel computing device is employed for achieving much higher statistics. Then 
one is able to look for and calculate transverse flow of each individual cell more accurately. 
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D. Results 

1. Rapidity distributions 

We now present first numerical results obtained for the time evolution of the gluons 
produced in a central Au+Au collision at RHIC energy ^/s = 200A GeV. In the simulations 
the number of the test particles is set to be Nt^st = 60. All results are obtained by an average 
over 30 independent realizations. Figures 1^ and 1^ show the particle number distributions 
per unit rapidity versus the space-time rapidity and the momentum rapidity at the times 
0.2, 0.5, 1.0, 2.0, 3.0, and 4.0 fm/c, respectively. The time interval of the overlapping 
for the two Au nuclei is about 0.17 fm/c. Therefore, the first extraction at 0.2 fm/c is 
just after the end of the production of the primary partons (or minijets). In Fig. I^one 
sees a noticeable spreading of the dN/drj distribution with progressing time. The reason is 
that the initially produced partons are distributed within a very small longitudinal region 
due to the Lorentz contraction of the Au nuclei. Their momentum rapidities, however, 
have a wider distribution, as can be seen in Fig. El The spreading of the space-time 
rapidity distribution continues until its width reaches a comparable magnitude with that 
of the momentum rapidity distribution. For the special case of a simple noninteracting 
free streaming system, the dN/drj distribution will then have exactly the same shape as the 
dN/dy distribution at late times. In the present case we see that at 4 fm/c the width of space- 
time rapidity distribution is about 4.2 and approaches nearly the width of the distribution 
of the momentum rapidity being about 5. It can be clearly seen that the spreading of the 
dN/drj distribution indeed slows down at late times. In the central space-time rapidity 
region the gluon density first decreases due to this spreading, and then increases because of 
the ongoing gluon production via the gg -^ ggg process. The gluon multiplication is most 
clearly demonstrated by inspecting the momentum rapidity distributions in Fig. |^ where 
for instance at y = the gluon number is double amplified until 4 fm/c. Moreover, at late 
times the net gluon production slows down, which implies the completion of the ongoing 
chemical equilibration. Of course, from the momentum rapidity distributions it is difficult 
to recognize any evidences for kinetic equilibrium. To investigate whether the system indeed 
does thermalize or not, one needs more detailed analyses in sufficiently local regions. We 
will present the results in next subsection. 
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Figure 021 shows the momentum rapidity distributions of the transverse energy (upper 
panel) and the total energy (lower panel) at the different times during the expansion. While 
the distributions would not change during an evolution like free streaming, we see in Fig. 021 
the decrease of the transverse energy and the energy transport from the center towards the 
higher rapidity due to the longitudinal work done by the pressure. This gives first significant 
indications of collective behavior. In addition, we note that when comparing Fig. ^2 with 
the upper panel of Fig. 021 the shape of the latter clearly looks more alike one dimensional 
Bjorken expansion than that for the particle number distribution. Hence, one cannot really 
conclude that the simple Bjorken expansion of constant dET/dy and dN/dy manifests. 

2. Thermalization 

In the following we study possible gluon thermalization in the "central region" being 
defined as a longitudinally expanding cylinder located in the middle of the expanding system. 
The radius of the cylinder is fixed to be 1.5 fm and its length is Arj = 1.0 from —0.5 to 0.5. 
In view of the possible buildup of transverse flow, one could consider a cylinder with varying 
radius which is comparable with the longitudinal length. On the other hand, however, 
the statistics within such cylinder would be very low at early times. Since the analysis of 
transverse flow, which we want to address in a further paper, shows that the transverse flow 
velocity is not large close to the central region even at time of 4 fm/c, the above choice with 
fixed radius is a reasonable compromise. 

In Fig. 021 we depict the time evolution of the gluon density and energy density in the 
central region. The densities are very high at early times. Alternatively, an implementation 
of the formation time for gluons, as briefly outlined at the end of Sec. IVIAI strongly reduces 
the densities at early times: ra ~ 20 fm~ and e ~ 50 GeV fm~ before 0.3 fm/c. After that 
time the results for the particle and energy density with and without the formation time are 
nearly identical throughout the subsequent evolution. It shows that even if the densities are 
very high at very early times with no formation time, the gluons rather stream freely within 
the first 0.3 — 0.4 fm/c. At 4 fm/c the energy density is still 1 GeV fm~^. Thus the parton 
picture of particle interactions is valid for the first 4 fm/c in a central Au+Au collision at 
RHIC After that hadronization should occur and the system is then in a parton-hadron 
"mixed phase" . 
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Figure E] shows the spectra of transverse momentum in the central region at different 
times during the expansion. The bold-folded histogram, which has a lower cutoff at 2 GeV, 
depicts the initial distribution of the primary gluons (minijets). The spectrum possesses a 
typical power-law behavior. Already at 0.2 to 0.5 fm/c a tremendous population of the soft 
gluons below 2 GeV has taken place. However, still a remedy of the edge at 2 GeV in spectra 
is visible. The "edge" vanishes at about 1 fm/c and the distributions become nearly exponen- 
tial and progressingly steepen at the later times 2, 3, and 4 fm/c. The ongoing steepening 
of the spectra in time represents a further strong indication of a (quasi)hydrodynamical 
expansion of an almost kinetically equilibrated system with decreasing temperature. 

To study kinetic equilibration in more detail, we first concentrate on the momentum 
anisotropy < p^ > /2 < p^ >. For an ideal, one dimensional boost-invariant hydrodynamical 
expansion the value of the anisotropy extracted within a region rj G [—Ari/2, A'q/2] is given 
by 



<Pt> 


pj_cosh(y-r)) 

Jl^dzJd^PTdyEple ^m 


2 <pl> 


Pi cosh(i/->7) 

2 Jl-dzld-^PTdyEple nr) 
/cf^/^dr^(coshr^)2/3 



6 Jo dr] (cosh r])^/^ — 5 /g di] (cosh r])"^/^ 
where z = t tanh(Ar7/2). The expression (ffn|) depends only on the longitudinal length of 
the local region where the momentum anisotropy is extracted, and goes to 1 in the limit 
Ar/ —>■ 0. In the central region with At] = 1, the anisotropy is equal to 0.65 for an ideal 
expansion. In Fig. |33 the time evolution of the momentum anisotropy extracted from the 
present simulations is depicted by the solid curve. Compared with the thermal value (0.65), 
the curve in Fig. 1^ shows first a significant increase during a short time 0.6 fm/c and then 
a smooth relaxation to that thermal value. The early increase of the momentum anisotropy 
is due to the initial pt cutoff, px > 2 GeV, and the fact that the primary gluons with large 
longitudinal momentum also have large rapidity and thus move rapidly out of the central 
space-time rapidity region. The following decrease of the anisotropy unambiguously implies 
the ongoing persistance of kinetic equilibration. The reason why the anisotropy is still larger 
than the thermal value is due to the fact that particles with larger pt equilibrate later, as also 
seen from the pr spectra in Fig. El From that particular analysis, quantitatively, the gluon 
system becomes approximately fully equilibrated at 2.5 fm/c. On the other hand, as just 
stated, the clear bending over at a time of 0.75 fm/c signals that the strong thermalization 
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has already started at that time, as one also notices from the onset of the pronounced 
exponential behavior at a similar time as seen in Fig. EJ 

The rapid streaming of the high-energy particles out of the central region at the beginning 
of the expansion also explains the dramatic decrease of the gluon density, energy density 
(both shown in Fig. I^Hj) and the effective temperature T = e/3n at early times, which 
is shown in the upper panel of Fig. H^l by the solid curve. The further decrease of the 
temperature until 200 MeV at 4 fm/c is due to the fact that work is done by the pressure 
and also due to the ongoing production of gluons. In case of simple free streaming the 
effective temperature would be constant over the whole time. To characterize the time 
dependence of the temperature we assume that the temperature behaves like T ~ 1/t" with 
a time dependent exponent. a{t) is shown in the lower panel of Fig. EHl We see that 
the exponent is almost constant, about 0.6, at late times and is roughly double the size 
of 1/3, which one expects for an ideal, one dimensional boost-invariant expansion. This 
might indicate the buildup of transverse flow, but is mainly due to the further production 
of gluons. For this we also extract the gluon fugacity from the simulations and depict its 
time evolution in Fig. El in a way similar as in Sec. IIVI [see Eq. ()38|l ]. Until 4 fm/c the 
chemical equilibration is still not fully achieved. Inspecting again Fig. 021 we have also 
plotted there for comparison the standard Bjorken behavior n ~ 1/t and e ~ 1/t^^^ with a 
fixed intercept at time t = 0.5 fm/c. One clearly recognizes that the particle number density 
decreases more slowly (with an exponent of about —0.7) due to the particle production. On 
the other hand, most interestingly, the energy density more or less exactly follows the form 
which one would expect from ideal Bjorken hydrodynamics. Indeed, the standard relation 
P = e/3 is all what enters into ideal hydrodynamical evolution for massless constituents, 
irrespective whether the system is chemically saturated or not. This finding gives another 
evidence that the system in the central region behaves nearly as an ideal fluid. We conclude 
that starting from the special, yet highly nonthermal initial condition a gluon plasma, even 
not fully thermalized, may form at 1 fm/c in a central Au+Au collision at RHIC energy and 
its ongoing evolution in bulk behaves (quasi)hydrodynamically. 

We have to note here that, of course, this reasoning will depend crucially on the initial 
conditions chosen. If we would only double the number of initial gluons, thermalization 
should roughly occur twice as fast. Indeed, our initial gluon number is lower compared to 

% ISQl , where a factor of 2 — 4 more initial gluons is assumed. 



other studies in the literature l2f 
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This will then clearly imply that full gluon equilibration within a consistent pQCD approach 
can have a full realization at RHIC. A detailed study, addressing various initial conditions 
for the gluon number, i.e., different forms of minijet productions or color glass condensate 
initial conditions, has to and will be done in a further publication. 

Figure EHl shows the time evolution of the cross sections which are first calculated as 
ensemble averages over all the possible collisions in a cell and then averaged over all the 
cells within the central region. As < Vrei >~ 1 in the central region, the collision rates of 
the gg *-^ gg and gg -^ ggg are obtained by -R = n < fre/O" >~ n < a >, respectively. 
We have compared these collision rates with those counted directly from the simulation and 
have seen nice agreements. The increase in time of the two cross sections is due to the 
fact that the cross sections are inversely proportional to the screening mass squared and 
the latter is proportional to the temperature squared. One sees that (Tgg^gg is always larger 
than cTgg^ggg- FoT kiuctic equilibration, however, not only a large total cross section but 
also large scattering angle are essential for a possible fast thermalization. In other word, the 
transport cross section 

at= da sin^ 9c.m. = / dOc.m. -7^ sin^ 9c.m. (71) 

is the key quantity controlling the ongoing of the equilibration by given particle density n. 
Oc.m. denotes the scattering angle in the center of mass frame of the colliding particles. For 
^ 99 ^^ 999 process each outgoing particle has its own scattering angle. In this case we 
modify Eq. (f7T|) by (sin^ 61 + sin^ 62 + sin^ ^3)/3 instead of sin^ Ocm.- The averaged transport 
cross sections are shown in Fig. UHlby the short dashed and short dotted curves. Taking 
into account that at late times the collision rate of the ggg —>■ gg is comparable with the rate 
of the gg — > ggg process, one realizes that the inelastic processes are actually the dominant 
processes driving the system to kinetic equilibrium. The ensemble averaged running coupling 
< a^ > is also extracted within the central region during the gluon evolution. It turns out 
that < a^ > increases almost logarithmically in time from 0.2 at 0.2 fm/c to 0.5 ~ 0.7 at 
4 fm/c. We note that when comparing the cross sections calculated in thermal equilibrium 
(see Fig. ITT|l . the cross sections Cgg^gg and agg^ggg extracted from the dynamical runs are 
4 ~ 5 times larger at later times. This is because first a^ had been fixed to 0.3 in Fig. ^2 
and is thus smaller. Second, the screening mass is appreciably smaller in the dynamical 
calculation as the gluons are not fully saturated in its occupation number. Both effects add 
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up to the difference. 

Following the expression of the differential cross section one knows that the gluon elastic 
collisions favor small angle scatterings. The transport cross sections in Fig. EHl indicate that 
the angular distribution of the inelastic collisions is more moderate than that of the elastic 
collisions. As can be realized from the differential cross sections expressed in Appendices O 
andini the angular distribution of the elastic scatterings depends on m|,/s, while it depends 
on m|,/s and A^^i for the inelastic collisions. In Fig. HHlwe depict the angular distributions 
of the gg -^ gg and gg -^ ggg scatterings for the parameters m'jj/s = 0.05 and Xg\^ = 4. 
The distributions are calculated according to the differential cross sections. The two param- 
eters are chosen from an intermediate situation within the simulation. We see that while the 
angular distribution of the elastic collisions clearly shows forward scatterings as expected, 
the angular distribution of the inelastic collisions is surprisingly almost isotropic. The reason 
for this behavior is due to the effective Landau- Pomeranchuk cutoff being implemented. For 
a larger Xg^/s the gg ^-^ ggg processes would also favor the more the small angle scatterings. 
Notice that 6*3 denotes the angle of the radiated gluon and thus possesses also a cutoff in 
its distribution due to the incorporation of the Landau-Pomeranchuk suppression of low 
momentum gluon emissions. 

In Fig. EOlwe present the px spectra at different times in the central momentum rapidity 
integrated now over the whole transverse region. The spectra are arranged in the same way 
as in Fig. EJ Comparing the spectra in Fig. EHlwith those in Fig. El we see that there is 
no full global thermalization over whole transverse region until 4 fm/c. At least for the lower 
momenta we see a nearly exponential population and a clear steepening at the later stages. 
Part of the minijet spectra, of course, survives as those gluons might escape directly from 
the outer region without interactions. In addition. Fig. EDI also demonstrates the potential 
energy loss of gluons due to the Bremsstrahlung process. The new developed parton cascade 
offers an alternative possibility to investigate the phenomenon of the jet quenching in a more 
quantitative way based on a full 3 + 1 dimensional treatment of the geometry. To be able to 
compare the numerical results with the experimental data one has to model the mechanism 
of the hadronization and include further hadronic interactions. A detailed analysis is again 
one of possible future projects. 
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3. Simulation only with elastic scatterings 

In order to further focus on the importance of the inelastic channels to the evolution, 
to the thermalization and to the potential onset of nearly ideal hydrodynamical behavior 
of the partonic system, we now perform simulations, for comparison, only with pure elastic 
scatterings among the gluons. Since in this case no gluons will be produced during the 
evolution, more test particles are needed to build for a fine cell structure. We set Ntest = 240. 
Fig. 1^ depicts the spectra of the transverse momentum in the central region at different 
times. The population of the soft gluons below 2 GeV is rather low and the distributions at 
large pt are only slightly altered. Indeed gluons with highest momenta get more populated. 
It is obvious that the gluon system is not thermalized during the expansion. This can also 
be realized from the time evolution of the momentum anisotropy presented in Fig. |331 by 
the dashed curve, where the anisotropy saturates at much higher value than 1 at late times. 
Furthermore, the constant temperature shown in Fig. E^l indicates that the evolution of 
the gluons is almost close to free streaming. (Please note that the abrupt decrease before 
0.5 fm/c is also due to the free streaming of the energetic gluons out of the central space- 
time rapidity region.) We remark that in the full dynamics with the inelastic channel, the 
contribution of the elastic scatterings to kinetic equilibration is actually significantly larger 
than that shown here, because in the full dynamics we have more gluons due to the radiation 
and the transport cross section also becomes larger compared to a nonequilibrated system. 

In principle, kinetic equilibration can be achieved by elastic scatterings alone, if ad hoc 
the transport cross section is chosen sufficiently large. To demonstrate this we carry out 
simulations with isotropic collisions and a large and constant total cross section of a22 = 30 
mb. The corresponding transport cross section is thus 20 mb. Such extreme conditions of an 
assumed large opacity in 2 <-> 2 reactions have been used in Ref. J28[ to study the possible 
buildup of elliptic flow. Figure EH shows the pt spectra in the central region at different 
times. Indeed we observe fast equilibration. The spectrum at 0.5 fm/c is almost thermal. 
At the later times the distributions become more and more steeper, which indicates the 
cooling down of the system due to (quasi)hydrodynamical expansion. The time evolution 
of the momentum anisotropy, the dotted curve in Fig. |^ shows that from 1.0 fm/c the 
anisotropy is almost constant and slightly higher than the value of the ideal, one dimensional 
boost-invariant expansion, 0.65. Moreover, also the exponent describing the cooling of the 
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temperature (see Fig. 1^]) . a{t), is nearly constant from 1 to 3 fm/c and only slightly 
greater than the value of the ideal expansion, 1/3. All this demonstrates that for the 
given extreme conditions the gluon system equilibrates indeed rapidly and then expands 
nearly hydrodynamically according to the ideal Bjorken scenario. However, of course, the 
constant and isotropic cross section cannot be further motivated. In addition, following that 
particular evolution, the system would stay for a rather long time in a hot, but very dilute 
and undersaturated (in its gluon number) deconfined state (see upper Fig. im|) . Contrary, 
in the more realistic situation with inelastic collisions included, the temperature drops much 
more dramatically and the system would stay only until t ^ 4 fm/c in a pure deconfined 
state, being then (nearly) fully saturated in the gluonic degrees of freedom, and will then 
hadronize. 

Figure I3H1 shows the time evolution of the transverse energy per unit momentum rapidity 
at midrapidity for the three cases compared also in Figs. 021 and EHl We see that the trans- 
verse energy decreases in the simulation including pQCD elastic and inelastic interactions 
and in the simulation employing an isotropic, large cross section. In contrast to the cooling 
of the temperature, to which the production of gluons also contributes, the decrease of the 
transverse energy within a unit rapidity is purely due to the longitudinal work done by the 
pressure! In the simulation employing large and constant cross section, energetic gluons 
are extremely stopped during their formation periode, so that the transverse energy is very 
large at very early times and pressure seems to be already built up during the overlap of two 
nuclei. This leads to the following (unrealistic) strong explosion with drastical cooling. The 
unaltered behavior of the transverse energy in the simulation including only pQCD elastic 
scatterings indicates again that in this case the parton evolution resembles a free streaming. 
One observes that from times t ^ 0.5 fm/c both the full pQCD (including gg ^^ ggg) and the 
"strongly coupled" (with isotropic (722 = 30 mb) evolution show almost the identical value 
and the same decrease in time for the total transverse energy per rapidity. This again mani- 
fests that both pathes resemble (quasi)hydrodynamical behavior by performing a significant 
amount of (longitudinal) work. 
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VII. SUMMARY AND OUTLOOK 

We have developed a new 3 + 1 dimensional relativistic transport model solving the kinetic 
on-shell Boltzmann equations. Besides binary 2 ^^ 2 scatterings, inelastic 2 ^^ 3 processes 
are also implemented in the cascade. The numerical emphasis is put on the extension of 
the stochastic collision algorithm for the back reaction 3 — i> 2 which is treated - for the 
first time - fully consistently within this scheme. Although the development specifically 
aims at a simulation of the parton evolution in an ultrarelativistic heavy ion collision, the 
presented algorithm will certainly have more potential applications beyond the scope of 
this paper. Also the standard geometrical collision algorithm (based on the geometrical 
intepretation of cross section) has been discussed in detail. In particular, we find out that 
for the case that the mean free path of particles is in the same order as or comes below the 
interaction length, which is always true in a very energetic (and dense) high-energy heavy 
ion collision, the results from the simulations employing the geometrical method have shown 
several unphysical numerical artifacts. The convergency of the numerical solution in this 
scheme for Ntest —>' oo turns out to be not as efficient as it does in the simulations when 
employing the stochastic method. 

The operation of the newly developed cascade has been demonstated by investigating 
gluon thermalization for a central Au+Au collision at RHIC energy. The numerical results 
show that starting initially from a nonthermal system made up of minijets (with cutoff 
Pt > Po = 2 GeV), the gluons in the expanding center equilibrate kinetically on a scale of 1 
fm/c and evolve further according to (quasi)hydrodynamics. The system cools down due to 
the hydrodynamical expansion and ongoing gluon multiplication. Full chemical equilibration 
follows on a longer time scale of about 3 fm/c. We have studied the contribution of the elastic 
and inelastic collisions to kinetic equilibration. It turns out that the inelastic scatterings 
are the main responsible processes driving the system to equilibrium. Without any inelastic 
collision channel, the collective behavior observed nowadays at RHIC cannot be generated, 
unless one uses an unrealistic large cross section (or equivalently a large gluon density) to 



mimic a strongly interactive gluon system 



3- 



We have also realized that the angular 



distribution of the gg ^-> ggg processes is almost isotropic during the expansion. This leads 
to larger transport cross section compared with the elastic scatterings. 

Even in the simulations applying the stochastic algorithm, particles do collide at nonzero 
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distance due to the nonzero spatial subvolume. Therefore, one may worry about acausal 



effects due to larger signal velocity than c in the cells |24l l37| . In principle, the spatial 
cell length should resolve the spatial inhomogeneities in the dynamical system. For the 
situation when using a 30 mb cross section for mimicking a strongly interacting system, 
the mean free path of the particles is initially smaller than the transverse cell length. To 
explore whether any potential acausal effect makes some numerical artifact, we now show 
a simulation employing half of the default transverse cell length and four times enhanced 
number of test particles (to keep the same statistics in cells). In Fig. |^we depict the time 
evolutions of the number and energy density of gluons extracted in the central region from 
the simulation with dx = dy = 0.25 fm and Ntest = 960 by the dotted lines. Comparing the 
results with those obtained with the default settings, depicted by the solid lines, we do not 
recognize any visible difference. This indicates that acausal effects seem to be not sensitive 
to the cell length when the system is rather uniformly distributed in space. 

Moreover, we note once more that the two timescales for kinetic and chemical equili- 
bration depend crucially on the initial gluon number. The one chosen here is rather low 
compared to other studies in the literature J2S|, |^, where a factor of 2 — 4 more initial 
gluons is assumed. This will clearly imply that the timescales for gluon equilibration within 
the present pQCD approach significantly shorten. Hence, In the future a lot of details have 
to be explored: Thermalization (also of the light and heavy quarks degrees of freedom) has 
to be investigated with various initial conditions like minijets, with a detailed comparison to 
data, or the color glass condensate, serving as input for the so-called bottom up scenario of 
thermalization. How likely is the latter picture for true coupling constants and not paramet- 
rically small ones? Furthermore, the indication of the hydro dynamical behavior during the 
expansion, which is one of the main findings from our first and exploratory study concerning 
RHIC physics, gives strong motivation for exploring transverse and elliptic flow using this 
new kinetic parton cascade. Can the inelastic interactions generate the seen elliptic flow f 2? 
Furthermore, one can also compare the present calculations with some fixed and specified 
hydrodynamical initial conditions directly with calculations based on viscous relativistic hy- 
drodynamics 59|, either assuming Bjorken boost invariance within an expanding tube or for 
full 3-1-1 dimensions. Such a comparison can tell how viscous the QGP really turns out to 
be. More practically, also the phenomenon of the jet quenching or electromagnetic radiation 
can be studied systematically within the new transport scheme. 
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Finally, the technique of the parallel programing is needed to improve the practical opera- 
tion of the cascade. With this technique quantities like the screening mass can be calculated 
and incorporated more precisely and quantum effects like the Pauli-blo eking and gluon en- 
hancement can be then implemented straightforwardly. 
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APPENDIX A: COLLISION TIMES IN THE GEOMETRICAL METHOD 

Within the algorithm implementing the geometrical picture collisions occur if the con- 
sidered particles approach each other and their closest distance is less than the interaction 



length Jojix. This criterion will be inspected in the center-of-mass frame of the colliding 
particles. Suppose that f, = (tj, Tj), -pi = {Ei, pj) and r^ = (t^, rQ, p'i = {E[, y>[), i = 1, 2, are 
the space-time coordinates and four momenta of two particles in the lab frame and in their 
cm. frame, respectively. Defining H = {r2—'f'i)-{pi+P2), one has in the cm. frame: t[ > t'2 if 
H < and t[ < t2ii H > 0. For the case t[ > t^ (otherwise we change the indices of the par- 
ticles) the two particles will approach each other if p^ [pi-{f2 — fi)\ — {pi-P2) [P2'(^2 — "^i)] < 0. 
The closest distance of the colliding particles in the cm. frame is 



where 



Ar:^f/- '"''X';/'"" - (^1' 



a = {f-2-ri)-pi, b= {f-2 - f-i) ■ p2 , 

c = p\, d = pl, e = pi-p2, 

f = if-2-hf. (A2) 
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If Ar'g < a/ct/tt, the particles will collide at the same time t'^i = t'^2 ^^ ^^e closest distance 

in the cm. frame. Making Lorentz transformation back to the lab frame gives 

ad — be bc — ae 

tci = ti — El — J , tc2 = ^2 + -c/2 -^ J • (A3) 

e^ — ca e^ — ca 

We call tci and tc2 the collision times. Due to the spatial separation, the two collision times 

have, in general, different values, tci 7^ ^c2- This means that one of the particles reacts 

later within the same collision. The transformed space coordinates at the collision times are 

correspondingly denoted by Td and rc2- The new momenta of the particles are sampled in 

the cm. frame according to the given differential cross section and then transformed to the 

lab frame, which are denoted by pd and Pc2- We thus label the particles with (td, rd) and 

{Eci,Pci), {i = 1,2), and keep the labels until their next respective collisions. For example, 

ti denotes the time when the last collision of particle 1 occurs. It is kinematically possible 

that the case ti < t^ < ^2 < tc2 occurs. Such a collision sequence is not causal, because 

at tci when the particle 1 experiences the collision with the particle 2, the particle 2 is just 

on the way to its last collision with some other particle. To forbid those collisions we add 

an additional criterion: The collision times td and tc2 should be greater than ti as well as 

^2- Illustratively, the additional criterion means that during the time interval |ti — ^2!, the 

particle, which will change its trajectory later (it is the particle 2 in the example), is not 

considered for dynamics for that particular interval. 

In the following we are interested in the probability distribution of the difference of 

collision times. Ate := \tci — tc2|, in a thermal system of massless particles. In this case we 

have c = d = 0. If ti 7^ t2 (e.g., ti < t2), the particle with smaller time (ti) can propagate 

freely to the larger time (t2), which does not give any effect on the whole evolution due to 

the additional criterion. Thus we obtain 

Ui +U2 



Ate = ri2 



(A4) 



1-u 

ri2 denotes |r2 — ri| and Ui = cosa^, u = cos 6*, where ttj is the angle between pj 
and r2 — ri and 6 is the angle between pi and p2. Since u relates Ui according to 



u = Jl — u{ Jl — U2 cos(0i — ^2) + Ml M25 where 0j is the polar angle of pj around r2 — ri, 
()A4j) can be now expressed by Ate = ''"12 -^(^i, ^2, 0) with := 0i — 02- One obtains the 
probability distribution of Ate by the integral 

P{Atc) = J_ duij_ du2J" d<pj dri2P(ri2,ui,M2,0)5(Ate-ri2F)e(y^- Ar^) 
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= / dui du2 d(j)P{ri2,Ui,U2,(p)\r^2=At,/F^r, tt ©(voT^ - Ar^) , (A5) 

j-i J-i Jo r [Ui,U2,(p) * 

where P{ri2,ui,U2,(p) is the muhiple probabihty distribution. Note that it is easy to re- 
ahze that Ar^ can also be expressed as a function of ri2, Ui, U2 and cf). Since ri2, Mi, 
U2 and are independent variables, P{ri2,Ui,U2,(f)) can be factorized, P{ri2,Ui,U2,4') = 
P{f'i2)P{ui)P{u2)P{(f)). For a thermal system we have P{ui) = 1/2, P(0) = l/27r and 
-P(^i2) = ^^12/P^i where R serves as a normalization factor and is set to be much larger 
than the interaction length. We realize that the probability distribution ()A5|) depends only 
on the size of the total cross section. For a constant cross section we calculate the integral 
in ()A5|1 numerically. Figure EHl shows the results for a = 10 mb and a = 30 mb. The 
distribution has a larger width for larger cross section. We also calculate the mean value of 
Ate and obtain < Ate >= 0.24 fm/c for cr = 10 mb and < Ate >= 0.41 fm/c for a = 30 mb. 

APPENDIX B: OPTIMIZATION OF THE COMPUTING TIME WITHIN THE 
GEOMETRICAL METHOD 

Consider a system with A^ particles in total. To get the next collision, A^(A^ — l)/2 
operations have to be carried out to get all the ordering times for each particle pair and 
N{N — l)/2 — 1 comparisons have to be made to obtain the particle pair which collides next. 
Then these two particular particles propagate freely until the two respective collision times 
when the respective momenta will be sampled according to the differential cross section. 
The same procedure will be repeated as long as needed. Since the operation number in 
each step is proportional to A^^, the computing time increases strongly with increasing 
particle number and increasing collision number. However, a large amount of operations are 
obviously futile, because after the update of two colliding particles only the ordering times 
of particle pairs which involve one of the two updated particles are indeed needed. Therefore 
only 2(A^ — 2), but not A^(A^ — l)/2, operations are necessary if one, in principle, wants to 
save all the ordering times from the last step. This, of course, reduce the computing time 
enormously. On the other hand, however, a large storage for those ordering times would be 
required. For an optimization we thus do not store all the ordering times, but only do store 
for each particle the informations of its possible next collision: the ordering time and the 
collision partner. We need therefore 2A^, instead of A^(A^ — l)/2, memory places. The next 
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collision will be obtained by comparing the marked and stored ordering times. In a next step 
we compute only the ordering times of the last colliding particles with the other particles 
{2{N — 2) operations) and compare them with the other stored times, respectively, to obtain 
the new informations of the next collision for each particle. The "worst" case then occurs 
when the next collision partner of a particle is one of the last colliding particles. In this case 
the stored informations for this particle are out of use and one has to compute the ordering 
times of the considered particle with all the other particles (additional A^ — 3 operations). 
Fortunately those cases do not happen frequently in practice. We note that our prescription 
is different from the optimization used by Zhang in his parton cascade j2^, which follows 
the fact that particles which are far away from each other most probably do not collide as 
next pair. In this algorithm the space was divided into cells and only particles from the 
same cell and the neighboring cells may collide next within the geometrical method. 

APPENDIX C: PARTON-PARTON SCATTERING CROSS SECTIONS 

Differential pQCD parton-parton cross sections in leading order of ag have been calculated 
in Ref. J48|. For elastic gluon scattering the differential cross section reads 

~^r~ ~ ~2^^'^ " ^ ~ ^ " ^^ ' ^^^^ 

where s, t and u are the Mandelstam variables, —t is equal to the momentum transfer 
squared 

-t = g2 = |(l-cos0), (C2) 

where 6 denotes the scattering angle in the cm. frame of colliding partons. For small angle 
scatterings the momentum transfer is approximately equal to its transverse component q±. 
Therefore we have —t ~ g^. Since the differential cross section ()C1|) diverges at small t 
(and also at small u due to the symmetry of identical particles), Eq. ()C1|) can be expressed 
approximately as 



da93^9g 9vra2 



(C3) 



dql {qlY ■ 

We regularize the infrared sigularity in Eq. ()(^3j) employing the Debye mass and obtain 

dql {ql + miy ^ ^ 
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dql {ql + miy 



The other approximate differential cross sections are achieved in the same way and read as 
follows: 

(C6) 
(C7) 
(C8) 
(C9) 






s 



dql %ql + mlY' 



dq]_ dq\ %(l\ + ^d)^ ' 

da^^^^g 647ra2 



dq\ 21s{q\ + m?) ' 



r,' Pi' 



dai^^q'q' ^-j^otl t^ + u^ 



where m|, and m? denote, respectively, the Debye mass for gluons and for quarks. In the 
last expression —t is not replaced by q\, since qq -^ q'q' processes do not favor small angle 
scatterings. Employing the fomulas ()C4|) - ()C10|) the total cross sections can be obtained 
analytically by integration. Equations ()C4|) - ()C10|) also then dictate how to sample new 
momenta for particles after an occurring collision. 

APPENDIX D: CROSS SECTION FOR gg ^ ggg PROCESSES 

For the multiplication process gg — > ggg^ the Gunion-Bertsch formula 49] is used for the 
matrix element squared in leading order of pQCD, and modified by implementing the Debye 
screening mass. This is 

_ /v s' \ I i2iA 1 (Di) 

where g^ = 47ras, and q^ and k_|_ are, respectively, the transverse component of the mo- 
mentum transfer and that of the momentum of radiated gluon in the cm. frame of the 
colliding gluons. In this section we will give the derivations of the cross section Cgg^ggg and 
Iggg^gg deflucd lu Scc. Illll bv au integral of the scattering amplitude given in Eq. ()D1|) over 
momentum space. 

Employing usual convention, the total cross section for a. gg ^ ggg process is defined as 

^"^ =2-8 J {2.f2E[{2.f2E',{2.Y2E'} ^''-"'"'\^^^''^ ^ \Pi + P2 ~ P. - P. - Ps) , 

(D2) 

58 



where pi, p2, p'l, P2, and p'^ are the four momenta and all momenta are expressed in the cm. 
frame of the two colliding gluons. We assume that 3' denotes the radiated gluon. Integrating 
over (Pp2 gives 

= ^^J <i'q^dy[d^k^dy\Mn^r2-3'\^ HF) , (D3) 

where y[ and y denote the momentum rapidity of 1' and 3', respectively, and 

F = {pi+p2-p'i -p'zf 

= s — 2y/sq± cosh y[ — 2\/sk± cosh y + 2q±k± cosh y[ cosh y 

+2q^ ■ k_|_ — 2g^/c^ sinhy^ sinhj/. (D4) 

Further integration over y[ gives 



^99^999 2567r^ 



\9y[ 



F=0 



where all the solutions of F = contribute to the sum. The corresponding differential cross 
section has the form 

"-^99^999 _ 1 I »^ |2 \^ -*- (-ptf;\ 



d'^q±d'^k±dy 2567r^s \^ 

\9yi 



F=0 

dF 



dy'i 



F=0 



This is different than in Ref. [3l|, where the authors ignored the factor J2^/ 
However, to make the correct implementation of the detailed balance for gg ^^ ggg processes, 
one should take the exact formula of the cross section. Expressing d'^q± and d'^k± in polar 
coordinates and integrating one of the two angles, one obtains 

/ d^q^d^k^dy -^ ii l dq\dk\dy / d(j) , (D7) 

where (f) denotes the angle between q_L and k^. 

We now turn to determine the integral boundaries for Eq. ()D7|) . At first, the energies 
of the three particles in the final state cannot be greater than \/s/2 because of the energy 
conservation. (Note that the total energy is equal to a/s.) We then have the upper bound- 
aries for q\ and k^: q\ < s/4 and kj_ < s/4. Secondly, k± and y will be further constrained 
by Q{k±Ag — coshy) due to the Laudau-Pomeranchuk suppression (compare Sec. II V|) . This 
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a~ / dql I dkl I dyl dcj)--- (DIO) 



leads to a lower cutoff for k±_: k± > 1/Ag. For given q± and k±, the constraints for coshy 

are 

E' fs 
cosli-u < k\ka and cosli-u = — ^ < — — . (D8) 

^ - ^ ^ k^- 2k^ ^ ' 

Thus the upper boundary of y, denoted by ?/„, is the smaller one among Arcoshik^Kg) and 
Arcosh{^\fsl2k^. Finally we have 

/•s/4 /•s/4 |-y™, /-TT 

^99-999 ~ / dq^\ dk^\ dy dcp--- . (D9) 

JO J1/A2 J-y^ Jo 

This integral actually scales with s, agg^ggg = a/s, where 

"1/4 ^ /.1/4 _^ py. 

-ym "'O 

with q\ = qj_/s, k^ = kj_/s, Ag = Agy/s, and m|, = m'jj/s. a depends on two parameters: 
Ag and rh'jj. We evaluate the above integral numerically using the Monte Carlo integration 
routine VEGAS [50|. For any sampled point {q\^k\^y^(f)) one has to solve y[ for F = 
in Eq. ()D4j) . If there is no solution, then the chosen point is out of the kinematic region 
and thus has no contribution to the integral. Thus the equation F = serves as a further 
constraint for the kinematic region of collisions. 

For the annihilation process ggg -^ gg, the analogous quantity as CTgg^ggg, which sums 
all the possible final states, is Iggg-^gg defined via 

•'«-» = \ I (afk (sitl' 1^— I^P^)'*"'"' + "^ + "- "'' - "» • (°"' 

where the factor 1/2 takes the identical gluons 1' and 2' into account and 

|A^123^1'2'P = — |A^l'2'^123r , (D12) 

where z/^, = 2 x 8 is the gluon degeneracy factor. Since Iggg-,gg is invariant under Lorentz 
transformations, we evaluate the integral in the rest frame of the three incoming particles. 
Therefore it is pi + P2 + P3 = {V^, 0). Integrating over d^p2 in Eq. (JD11|) we find 

1 r d u' 

hgg-^99 = Y^J^\Ml23^l'2'fS{{pi+P2+P3-p[f) 

= -^ [ dE[d COS 9 d(l) E[\Mi23-.i'2'\^ 5(3 - 2./:s E[) 

= 7^ /' ^cos^ r d(P\Mi23^r2'\\ (D13) 

o47r^ J-1 Jo 
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where the sohd angle (cos^,0) defines the orientation of p'^ (p2 = — Pi because of the 
energy-momentum conservation). We now express q±, k± and q_L ■ k_L in |7Vfi23^i'2'P with 
the sohd angle and momenta of the incoming gluons. To do that one has to specify the 
fussion process 2^1. There are in total 6 combinations. Each combination contributes to 

Iggg^gg- Oue of them Is 

123 -^ 1'2' = (a) 23 -^ 2* and (b) 12* -^ 1'2' 

and p3 corresponds to (k_L,t/). In this particular case one can establish a coordinate system 
in momentum space whose Z axis coincides with the orientation of pi. We find out (after a 
direct but lengthy calculation) 

q^ = Eisme, (D14) 



k^ = E3JI — (sin 7 sin^ cos0 + cos 7 cos^)^ , (D15) 

q^ ■ k_|_ = -E1-E3 (sin 7 sin 6* cos^ cos0 — cos 7 sin^6') , (D16) 

where 7 denotes the angle between pi and p3, and (cos^,0) is, as defined before, the 
solid angle of p'^. Due to the Laudau-Pomeranchuk suppression Q{k±Ag — coshy) and 
coshy = E^/k± we obtain the kinematic region for the ggg -^ gg process 

k, > ^^ (D17) 

In analogy to agg^ggg, Iggg^gg sAso scalcs with s, Iggg^gg ~ l/s, wherc / depends on five 
parameters, namely Ei/ ^/s, E^/ ^/s^ 7, A^y^, and vn\)l s. 

APPENDIX E: MONTE CARLO SAMPLING OF MOMENTA FOR OUTGOING 
PARTICLES 

Momenta of outgoing particles are sampled in the rest frame of the incoming particles. 
Their momentum in the lab frame is obtained by Lorentz transformations. 

1. 2 ^^ 2 processes 

In the rest frame the energy of each particle is ^/s/2. The only to be sampled quantity 
is the solid angle (cos^,0). The scattering angle Q is samlped according to the differential 
cross section and the polar angle is sampled uniformly within [0, 27r]. 
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Since the pQCD differential cross sections ()C4|) - ()C10|1 can be integrated analytically, we 
can perform samplings for g^ (or cos 6') using the "transformation method" J50| from a 
uniform probability distribution. For isotropic collisions we sample the scattering angle 9 
according to the uniform distribution of cos 9. 

2. gg '^ ggg processes 

As shown in Appendix El the differential cross section for a gg ^ ggg process has the 
form 



E-^^^ (El) 



ay[ 



F=0 



dq\dk\dyd(j) {q\ + in\)Y ^^[(k^ — q±)^ + "^d] 

where (j) denotes the angle between q_L and kj_. We then first sample g_L, A;_l, y, and (j) 
according to Eq. (jElj) . Since the differential cross section cannot be integrated analytically, 
one cannot make samplings by means of the transformation method, as done for 2 — >■ 2 
processes. Instead, we employ the "rejection method" [50|. 

To make enough efficient samplings, we want to find out a special function of g^, /c^, 
y and 0, which should be always greater than the right hand side of Eq. ()E1|) at every 
point set (g_L, /c_l, y, 0) in the kinematic region and, more important, can be integrated out 
analytically over g_L, k^, y, and (p. Such a function is called as a "comparison function". If 
one has the comparison function, one can first use the transformation method to generate the 
random numbers according to the comparison function. Then one needs a further uniform 
sampling between zero and the value of the comparison function at the particular sampled 
point. If this random number is less than the value of the real distribution [right hand side of 
Eq. (jElj) ] at the sampled point, then we accept this sampling, if not, we reject this sampling 
and start a next trial. One possible choice of the comparison function is 

111 /^ X 

-9 9-79-^9""^) (E2) 

qj_ + mf) kj_ mjy 



dF 



dy[ 



F=0 



where m denotes a constant with a sufficient large value, which is greater than J2 1/ 

in Eq. ()E1|) at every point in the possible kinematic region. Since, unfortunately, one cannot 

obtain the upper limit for I] 1/ f4- , the value of m is an empirical number. 

We have to note that for an individual sampling one has to solve the equation F = 0, 
Eq. ()D4j) . Therefore one also obtains y[, the momentum rapidity of the particle 1', at the 
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same time when g^, k±, y, and are sampled. One sampling remains: The direction of 
qx is sampled uniformly in the transverse plan being perpendicular to the scattering axis. 
Finally we obtain the momenta of the outgoing particles 

Pi± = -q± > P'lz = (1± sinh y[ , (E3) 

P3± = k± , pL = k± sinh y , (E4) 

p'2 = -(p'i + P3). (E5) 

For a ggg -^ gg process the solid angle (cos 6, 0) is sampled again by using the rejection 
method. We find 



^■^ggg^gg ^ 1± 



dcosOdcj) (gi + "^i))^ fci[(ki — q±)^ + "^d] 

1 A„ 1 1 A„ 1 

ql + ml E3 ml Ef{l - cos2 6) +mlE^ml' ^ ' 

where we have employed the constraint ()D17|) and the identity ()D14|) for the particular 
example presented in Appendix iDl 



3. Isotropic 2 ^^ 3 processes 

A 2 ^ 3 process, 12 -^ 1'2'3', is assumed to be composed of a two-body scattering 
12 -^ 1'2* and a decay 2* -^ 2'3', where 2* denotes an intermediate state with an invarian; 



mass of m* = J El — p^. We employ the formula for the phase space integrations of 
and obtain the differential cross section of an isotropic collision 

ddO'i , 1 , ^^ /' ,„/ Er, 



where fii = (cos6'i,0i) denotes the solid angle of p'^ with respect to the collision axis, and 
Q2 = (cos 6^2, ^2) denotes the sohd angle of P2 with respect to p*, and 



f{E'^) = E,-E'^- ^pI + E'l - 2p,E!, cos 62 , (E8) 

A(s, ml, ml) = s^ - 2s{ml + ml) + {mj - ml)"^ . (E9) 



The integral over E2 in Eq. ()E7|) gives 

da23 {s - m1)ml 



dVLidm'ldVL2 {E^ — p* cos 62)' 
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(ElO) 



Qi, ml and Q2 are the to be sampled quantities. From Eq. ()E10|) we realize that the 
differential cross section does not depend on Qi and 02- Thus they are sampled uniformly. 
Integral over Qi and Q2 gives the probability distribution for ml. It is simply proportional 
to s — ml. We sample m^l by employing the transformation method. For given Qi and m^ 
the momenta of 1' and 2* are fully determined due to the energy-momentum conservation. 
Now Eq. ()E10|) just represents the probability distribution of cos6'2 for a given Qi and m,l. 
Its numerical sampling is straightforward. 

Sampling for an isotropic 3 — i> 2 process is more trivial, since just one solid angle is to be 
sampled and its probability distribution is uniform. 
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1 ' 1 ' 

Box calculation, L=5 fm, N=2000 
isotropic collisions, 0^2=10 mb 

geometrical method, N^^^, =1 , 

averaged over 50 runs 

analytical, T=2 GeV 




FIG. 1: Energy distribution at final time (t = 5 fm/c) of a system consisting of A'^ = 2000 massless 
particles in a fixed box. The initial energy distribution is set to be a delta-function at 6 GeV. 
The size of the box is 5 fm x 5 fm x 5 fm. We here apply the geometrical collision algorithm. 
The collisions are taken as isotropic and the total cross section is fixed to be a constant CJ22 = 10 
mb. The dotted line denotes the analytical result of temperature T = 2 GeV. The numerical 
distribution is obtained from an ensemble of 50 independent realizations. 
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FIG. 2: Collision rates for given particle densities. The size of the box is 5 fm x 5 fm x 5 fm. 
We apply here the geometrical collision algorithm. The collisions are isotropic and the total cross 
section is fixed to a constant (T22 = 10 mb. The particle system is taken initially as thermal with 
a temperature of T = 1 GeV. The solid line shows the expected relationship between collision rate 
and particle density: R = na22- The solid squares show the calculated collision rates without test 
particles {Ntest = 1) and the open squares show the results with 50 test particles per real particle 
(Ntest = 50). 
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FIG. 3: Energy distributions from box calculations. The thin histogram shows the same distri- 
bution as in Fig. ^ The thick histogram shows the result with a smaller total cross section of 
(T22 = 0.1 mb. 
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FIG. 4: Time evolution of the momentum anisotropy from box calculations. The initial condition 
and parameters are set to be the same as in Fig. ^ The dotted (dashed) curve shows the results 
obtained by employing the geometrical method without test particles (with 50 test particles per 
real particle). The solid curve shows the result obtained by employing the stochastic method with 
ten test particles per real particle. 
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FIG. 5: Energy distribution from box calculations. The initial conditions and parameters are set 
to be the same as in Fig. ^ We apply here the stochastic collision algorithm. The box is divided 
into equal cells. The length of a cell is 1 fm. 



71 



25 



20 



S 

i_ 

c 
o 

"o 
o 



— ' 1 ' 1 ' r 

Box calculation, L=5 fm 
isotropic collisions, cj,,=10 mb 

stochastic method 
■ N, =1 

test 




n[fm"1 

FIG. 6: Collision rates for given particle densities. The initial conditions and parameters are set to 
be the same as in Fig. |2I We apply here the stochastic collision algorithm. The cell configuration 
is the same as in Fig. |SJ 
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FIG. 7: Time evolution of momentum anisotropy from box calculations. The initial conditions and 
parameters are set to be the same as in Fig. ^ (or Fig. EJ. The stochastic method is used here. 
The cell configuration is the same as in Fig. [21 The curves show the results with different number 
of test particles. 
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FIG. 8: Time evolution of the particle density from box calculations. The initial conditions and 
parameters are set to be the same as in Fig. ^(or Fig. EJ. We consider isotropic inelastic collisions 
(2 <-^ 3) with a constant cross section of (T23 = 10 mb and employ the stochastic collision algorithm. 
The cell configuration is the same as in Fig. El The dotted line denotes the estimate using a simple 
time relaxation approximation. 
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FIG. 9: Energy distribution from the same calculations as in Fig. |H1 The histogram shows the 
numerical result. The dotted line shows the analytical expectation and the dashed line shows the 
analytical distribution (the same as in Fig. [S)) if the particle number would be conserved. 
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FIG. 10: Time evolution of the fugacity [n{t)]/neq versus the momentum anisotropy from the same 
calculation as in Fig. |B1 
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FIG. 11: Gluon collision rates and thermally averaged < fre/f > as function of temperature. 
The solid, dashed, dotted and dash-dotted line show the temperature dependence for gg -^ ggg, 
gg -^ gg, gq -^ gq and gg -^ qq transitions respectively. We consider here two quark flavors and 
employ a constant coupling a^ = 0.3 (for the cross sections and the screening masses). 
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FIG. 12: Time evolution of the gluon and quark number in box calculations. We consider here 
gluons and quarks with two flavors as parton species. Collision processes are the elementary two- 
body parton-parton scatterings and three-body processes gg ^^ ggg in leading order of perturbative 
QCD. The coupling is assumed to be a constant of Og = 0.3. The initial momentum distribution 
of particles is taken from the minijets production in central rapidity interval y £ (—0.5 : 0.5) in 
a nucleon-nucleon collision at RHIC energy y/s = 200 GeV. The initial particles are gluons and 
distributed homogenously in the box. The size of the box is 3 fm x 3 fm x 3 fm and the box is 
divided into equal cells. The length of a cell is 1 fm. The initial gluon number is set to be 500. 
The results are obtained from an average over 60 runs. 
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FIG. 13: Energy distributions at different times from the same calculation as in Fig. 1121 
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FIG. 14: Time evolution of the momentum anisotropy for gluons and quarks from the same calcu- 
lation as in Fig. 1121 
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FIG. 15: Time evolution of the temperature for gluons and quarks from the same calculation as in 
Fig. II2I 
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FIG. 16: Time evolution of the fugacity for gluons and quarks from the same calculation as in Fig. 

Ha 
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FIG. 17: Time evolution of the screening mass for gluons and quarks from the same calculation as 
in Fig. ini 
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FIG. 18: Time evolution of the collision rates for gluons in the different channels. The results are 
obtained from the same calculation as in Fig. 1121 
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FIG. 19: One dimensional expansion in a tube. The lab frame is labeled by X, Y, and Z, the 
boosted frame by X' , Y' , and Z' which is moving with a velocity of /? relative to the lab frame. 



85 



2000 




FIG. 20: Space-time rapidity distributions at different times {t = 0.11, 0.13, 0.16, and 0.2 fm/c 
from histogram with smallest amplitude to histogram with largest amplitude) from a simulation of 
one dimensional expansion in a tube. We consider a thermal and boost-invariant initial condition 
for evolving particles: Particles are produced initially on a hyperbola of tq = 0.1 fm/c and are 
distributed homogenously within a space-time rapidity interval rj G [—3 : 3], dN/drj^To) = 1748, 
which is depicted by the dotted straight line. The initial temperature is set to be T{to) = 2.6 GeV. 
The radius of the tube is i? = 5 fm. We consider 2 ^^ 2 collisions with isotropic cross section and 
a constant total cross section of CJ22 = 10 mb. The stochastic method is used in the simulation. 
The Tj bins of the cell configuration are set to be Arjc = 0.2. No test particles {Nfest = 1) are used 
in the simulation. The distributions are obtained by an average over 10^ independent realizations. 
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FIG. 21: Space-time rapidity distributions at time t = 0.2 fm/c in tube calculations. The initial 
condition and collision cross section are the same as in Fig. 1201 The stochastic method is employed 
in the simulations. The result showing structure with larger(smaller) period is obtained from the 
simulation with Aijc = 0.2(0.1). In the simulation with Arjc = 0.1 we use 2 test particles per real 
particle in order to achieve the same statistics in each cell as that in the simulation with Aijc = 0.2 
and Nfest = 1- The histogram, which is nearly constant, is obtained from the simulation with 
improved moving cell configuration of Ar/c = 0.2 and Nfest = 1- The dotted line shows the initial 
distribution dN/drj = 1748. All the distributions are received by an average over 10^ independent 
realizations. 
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FIG. 22: Time evolution of the particle density, energy density, and temperature extracted in the 
central space-time rapidity region r] £ [—0.5 : 0.5] from simulations of one dimensional expansion 
in the lab and boosted frame of a tube. The geometrical method is employed in the simulations. 
The initial condition and collision cross section are the same as in Fig. I2UI No test particles 
[Ntest = 1) are used. Only 2 <-^ 2 processes are included. The results are obtained by an anerage 
over 20 independent realizations. The thin lines indicate time evolutions of the quantities in the 
hydrodynamical limit. 
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FIG. 23: Time evolution of the particle density, energy density, and temperature extracted in 
the central space-time rapidity region rj £ [—0.5 : 0.5] from simulations employing the stochastic 
method in the lab and boosted frame of a tube. The initial condition and collision cross section 
are the same as in Fig. 1201 No test particles {Ntest = 1) are used. Only 2 •s-^ 2 processes are 
implemented. We apply the moving cell configuration with A?7c = 0.2. The results are obtained 
by an average over 20 independent realizations. 
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FIG. 24: Collision rate in the central space-time rapidity region for various particle densities 
experienced during the expansion. The results are extracted from the same simulations performed 
for the extractions of n(t) and e(t) in Figs. |221 and 1231 The solid line shows the analytical 
expectation. 
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FIG. 25: Averaged difference in space-time rapidity of colliding particles extracted in the central 
space-time rapidity region for various particle densities experienced during the expansion. The 
results are extracted from the same simulations performed for the extractions of n(i) and e(t) in 
Figs. 1221 and l23l The labeling of the symbols is identical to that of Fig. 1241 
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FIG. 26: Particle distributions versus space-time rapidity at the proper time r = 0.2 and 1.0 
fm/c extracted from simulations employing the geometrical and stochastic method in the lab and 
boosted frame. The initial condition and collision cross section (and cell configuration) are the 
same as in Fig. |22l (and in Fig. I2,3|) . In order to compare the distributions in the same physical 
regions directly, we have shifted the distributions in the boosted frame by — t/q = —2. Except 
that the distributions extracted from the simulations in the boosted frame using the stochastic 
method are obtained by an average over ten independent realizations, all other distributions are 
obtained from 20 independent realizations. The thin solid lines indicate the initial distribution 
dN/dr]{To) = 1748. The cut at r/ = 4 in the distributions at r = 1.0 fm/c for the expansion in 
the boosted frame is due to the fact that the end time of the simulation in the boosted frame is 
t' = 210 fm/c and thus particles with rj being greater than 6 (or 4 after the shift) have smaller 
proper time than 1 fm/c. 
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FIG. 27: Particle distributions versus momentum rapidity at the proper time r = 0.2 and 1.0 
fm/c extracted from simulations employing the geometrical and stochastic method in the lab and 
boosted frame. The results are obtained from the same simulations performed for the extractions 
of dN/dr}{T) in Fig. 1261 The thin solid curves indicate the initial distribution at tq = 0.1 fm/c. 
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FIG. 28: Distributions of the transverse momentum per unit rapidity at y = at r = 1.0 and 4.0 
fm/c (from upper to lower histogram) in a simulation employing the geometrical method in the 
lab frame. The initial condition and collision cross section are the same as in Fig. 1221 The results 
are obtained by an average over 20 independent realizations. The solid lines show the analytical 
distributions (|5H|) with the temperatures read off from Fig. 1221 
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FIG. 29: Distributions of the transverse momentum per unit rapidity at y = at r = 0.2, 1.0, and 
4.0 fm/c (from upper to lowest histogram) in a simulation employing the stochastic method in the 
lab frame. The initial condition, collision cross section, and cell configuration are the same as in 
Fig. 1231 The results are obtained by an average over 20 independent realizations. The solid lines 
show the analytical distributions (|5H|) with the temperatures read off from Fig. I23L 
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FIG. 30: Proper time evolution of the transverse energy per unit momentum rapidity at y = in 
the simulations employing the geometrical and stochastic method in the lab and boosted frame. 
The initial condition and collision cross section (and cell configuration) are the same as in Fig. [211 
(and in Fig. I23|) . The results are obtained by an average over 20 independent realizations. The 
thin solid line shows the analytical evolution in the hydrodynamical limit. 
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FIG. 31: Relative frame dependence of the particle density, energy density, and temperature in the 
simulations employing the geometrical method. The initial condition and collision cross section 
are the same as in Fig. 1221 The results are obtained by averaging 20, 2, and 20 independent 
realizations for increasing test particles Ntest = l, 4, and 25, respectively. 
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FIG. 32: Comparison of the space-time rapidity distribution with Ntest = 25 with the distribution 
without test particles at r = 0.2 fm/c. The distributions are extracted from the simulations 
employing the geometrical method in the lab frame by averaging 20 independent realizations. The 
initial condition and the collision cross section are the same as in Fig. 1221 The thin solid line 
indicates the initial distribution dN/dr]{TQ) = 1748 within ij £ [—3 : 3]. 
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FIG. 33: Convergency of temperature in the simulations in the lab frame with increasing test 
particles. The initial condition and collision cross section (and cell configuration) are the same as 
in Fig. |^(and in Fig. I23() . The results in the simulations employing the geometrical method are 
obtained by averaging 20, 2, 20, 5, 5, and 5 independent realizations for Ntest = 1, 4, 25, 100, 400, 
and 900, respectively. The results in the simulations employing the stochastic method are obtained 
by averaging 20, 10, 2, and 1 independent realizations for Ntest = 1; 2, 4, and 8, respectively. 
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FIG. 34: Collision rate and averaged difference in space-time rapidity of colliding particles. The 
results are extracted in the central space-time rapidity region rj G [—0.5 : 0.5] for various particle 
densities experienced during the expansion. The simulations are the same as performed in the 
upper panel of Fig. 1331 when discussing the convergency of the temperature. 
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FIG. 35: Time evolution of the particle density, energy density, and temperature extracted in the 
central space-time rapidity region rj £ [—0.5 : 0.5] from the simulations employing the stochastic 
method in the lab and boosted frame. The initial condition and cell configuration are the same as 
in Fig. [2S1 No test particles {Ntest = 1) are used. 2 •s-^ 2 as well as 2 <-^ 3 processes are included 
in the simulations. We consider isotropic collisions with constant cross sections of (T22 = 5 mb and 
(723 = 2.5 mb. The results are obtained by an average over ten independent realizations. The thin 
solid lines indicate time evolutions in the hydrodynamical limit. 
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FIG. 36: Time evolution of fugacity extracted from the same simulations performed for the extrac- 
tion of n{t) and e(t) in Fig. 1351 
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FIG. 37: Collision rate in the central region for various particle densities experienced during the 
expansion. The results are extracted from the same simulations performed for the extractions of 
n{t) and e(t) in the lab frame in Fig. |2S1 The solid squares, solid circles and open circles depict, 
respectively, the collision rates for 2 ^^ 2, 2 ^ 3, and 3 — > 2 transitions. The solid and dotted line 
show the analytical expectations. 
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FIG. 38: Particle distributions versus space-time rapidity and momentum rapidity at the proper 
time T = 0.2 and 1.0 fm/c, extracted from the simulations employing the stochastic method in 
the lab and boosted frame. The initial condition, collision cross section, and cell configuration are 
the same as in Fig. EH The distributions extracted in the lab (boosted) frame are obtained by 
averaging 20(6) independent realizations. The thin solid lines indicate the initial distributions at 
To = 0.1 fm/c. 
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FIG. 39: Convergence of energy density in the simulations in the lab frame with increasing test 
particles. The cascade simulations are performed employing the stochastic algorithm. The dotted 
lines depict the results with o"22 = 5 mb and o"23 = 2.5 mb, while the thin solid lines depict the 
results with purely elastic collisions and (T22 = 10 mb. The results are obtained by averaging 20, 
10, 2, and 1 independent realizations for Nfest = Ij 2, 4, and 8, respectively. The thin dashed lines 
show the hydrodynamical limit. 
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FIG. 40: Gluon number distribution versus space-time rapidity at the time t = 0.2, 0.5, 1.0, 2.0, 
3.0, and 4.0 fm/c during the expansion in a real, fully 3D central Au+Au collision at the maximal 
RHIC energy. 
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FIG. 41: Gluon number distribution versus momentum rapidity at the time t = 0.2, 0.5, 1.0, 2.0, 
3.0, and 4.0 fm/c during the expansion. 
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FIG. 42: Momentum rapidity distributions of the transverse energy (upper panel) and the total 
energy (lower panel) of gluons at the time t = 0.2, 0.5, 1.0, 2.0, 3.0, and 4.0 fm/c during the 
expansion. 
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FIG. 43: Time evolution of the gluon density and energy density in the central region: radial 
transverse extension xy < 1-5 fm and t] G [—0.5 : 0.5] for a central Au+Au collision at the maximal 
RHIC energy. The dotted curves denote the ideal hydrodynamical limit with a fixed intercept at 
time t = 0.5 fm/c. 
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FIG. 44: Transverse momentum spectrum in the central region at different times (t = 0.2, 0.5, 1, 
2, 3, and 4 fm/c from second upper to lowest histogram) during the expansion. The most-upper 
and bold-folded histogram with a lower cutoff at pT = 2 GeV denotes the spectrum of the primary 
gluons (minijets). 
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FIG. 45: Time evolution of the momentum anisotropy extracted in the central region. The solid 
curve shows the result from the simulation with full dynamics, while the dashed curve shows the 
result from the simulation with only elastic scatterings. The dotted curve depicts the result from 
the simulation with isotropic elastic collisions and with (unrealistic) large cross section of cj = 30 
mb. 
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FIG. 46: Time evolution of the effective temperature (upper panel) and the exponent describing 
the cooling of the system (lower panel) in the central region. The curves are arranged in the same 
way as in Fig. 1451 
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FIG. 47: Time evolution of the gluon fugacity extracted in the central region. 
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FIG. 48: Time evolution of the averaged cross section and the averaged transport cross section in 
the central region. The solid and dashed (short-dashed and short-dotted) curves depict the averaged 
cross sections (transport cross sections) for the gg -^ gg and gg -^ ggg processes, respectively. 
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FIG. 49: Angular distribution of the scattering processes gg -^ gg (solid curve) and gg -^ ggg 
for a representative situation during the gluon evolution. ^3 denotes the scattering angle of the 
radiated gluon and its radiation partner has the angle ^2- The distributions are computed with 
the parameters m?£)/s = 0.05 and Xg^/s = 4 extracted in the central region at an intermediate time 
during the evolution. 
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FIG. 50: Transverse momentum spectrum in the central space-time rapidity slice (?? G [—0.5 : 0.5] 
and all gluons in the transverse plan are counted for) at different times (t = 0.2, 0.5, 1.0, 2.0, 3.0, 
and 4.0 fm/c from second upper to lowest histogram). The most-upper and bold-folded histogram 
with a lower cutoff at px = 2 GeV denotes the spectrum of the primary gluons (minijets). 



116 



^ 10-^ 



1 ' r 







pQCD 

gg <-> gg 




3 4 5 6 

Pt [GeV] 



8 



FIG. 51: Transverse momentum spectrum in the central region extracted from the simulation 
with only elastic collisions at different times. The most-upper and bold-folded histogram with a 
lower cutoff at pT = 2 GeV denotes the spectrum of the primary gluons (minijets). According to 
the increase of the population of the soft gluons below 2 GeV, the other histrograms present the 
spectrum at times 0.2, 0.5, 1.0, 2.0, 3.0, and 4.0 fm/c, respectively. 
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FIG. 52: Transverse momentum spectrum in the central region extracted from the simulation with 
isotropic elastic scatterings and a large cross section of o" = 30 mb at different times (t = 0.2, 0.5, 
1.0, 2.0, 3.0, and 4.0 fm/c from second upper to lowest histogram). The most-upper and bold- 
folded histogram with a lower cutoff at pT = 2 GeV denotes the spectrum of the primary gluons 
(minijets). 
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FIG. 53: Time evolution of the transverse energy per unit momentum rapidity at midrapidity. The 
curves are arranged in the same way as in Fig. 1451 
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FIG. 54: Time evolution of the number (left panel) and energy (right panel) density extracted in 
the central region from the simulation with dx = dy = 0.25 fm and Ntest = 960 by the dotted lines, 
compared with the results with the default settings dx = dy = 0.5 fm and Attest = 240, depicted 
by the solid lines. 
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FIG. 55: Probability distribution of difference in "collision times" within the geometrical collision 
algorithm. In the calculations a thermal system is assumed and the cross section is set to be a 
constant. 
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